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Time series models with autoregressive, moving average 
and mixed autoregressive-moving average correlation struc- 
ture and with positive-valued non-normal marginal distribu- 
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leads to the formation of a likelihood function and a 
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Second, three methods for generating first-order 
moving average sequences with Exponential marginals are 
examined. These generalize the EMA(1) Exponential model. 
Negative correlation using antithetic variables is investi- 
gated in the moving average models. 

A preliminary analysis of wind speed data obtained over 

a 15 year period in the Gulf of Alaska is presented. A 
model with four harmonic deterministic mean multiplying 
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i IN@RODUCTION 


The classical approach to time series analysis based on 
linear, additive models with normally distributed, constant 
variance residuals is probably best presented by the work of 
Box and Jenkins [Ref. 1]. Although their work is widely ac- 
cepted and used, it is not applicable to some important time 
series. This is mainly because the Box-Jenkins approach is 
based on an assumed normal distribution for the series in 
question. However, the assumption of normality is not appro- 
priate when the series is known to be non-negative. Such 
series typically involve times between successive events in 
event processes. Examples are easy to construct. Times be- 
tween arrivals at a hospital emergency room, times between 
breakdowns in a tank main drive assembly, and times between 
detections of enemy armor vehicles are a sample of series of 
this type. Because of the non-negative nature of the series, 
the Box-Jenkins distributional assumptions and, hence, the 
analysis techniques are inappropriate. There is, of course, 
the possibility of data transformations but this is not appro- 
priate with very skewed marginal distributions and it is, in 
most cases, difficult to ascertain what the transformation 
does to the correlation structure of the series. 

Gaver and Lewis [Ref. 2] wrote the pioneering paper on 
the subject of autoregressive processes with non-normal marginal 
distributions. They presented the method for determining the 
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distribution of the innovative terms in the basic, linear, 


additive, autoregressive equations (first-order stochastic 


different equation) 


K = oX + € ( Lee) 


that was required to produce a given marginal distribution 
for the cosa sequence. They presented results for or se- 
quences with Exponential, Gamma, and mixed Exponential 
marginals. They also showed that this problem was the same 
as that of determining the class of self-decomposable (Type 
L) random variables (Feller, [Ref. 3], Loeve [Ref. 4]) al- 
though the connection between the solution to the self-decom- 
posable problem and equation (I.1) was not explicit in the 
literature. 

The Gaver and Lewis paper was followed by other papers 
which extended these results. Lawrance and Lewis [Ref. 5] 
presented a first-order moving average process with Exponential 
marginals. Jacobs and Lewis [Ref. 6] propounded a mixed auto- 
» regressive-moving average of order one, EARMA(1,1), and 
Exponential marginals. This was extended to an arbitrary 
order EARMA(p,q) process by Lawrance and Lewis [Ref. 7]. A 
further refinement of the first-order, Exponential, auto- 
regressive process (NEAR(1)) was presented by Lawrance and 
Lewis [Ref. 8]. While this contained the previous EAR(1) 
model, it did not suffer from the degeneracy inherent in (I.l). 


Jacbos applied these models to closed cyclic queueing networks 
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(Ref. 9] and Lewis and Shedler applied them to models of 
computer processes [Ref. 10]. 

This paper extends the results of these researchers and 
others in three areas. In Chapter II a mixed ARMA(p,q) model 
with Gamma marginals proposed by Lewis [Ref. 11], the GLARMA(p,q) 
model, 1S examined. The correlation structure is derived for 
peweral values of p and q. Of particular note is the AR(1) 
case (p = l, q = 0), called GLAR(1), where the conditional 


density of x given X is derived. This leads to the deri- 


n-1 
vation of a likelihood function and a numerical technique to 
evaluate and maximize the likelihood function with respect to 
the model parameters. This provides a useful technique for 
estimating model parameters. Using this numerical technique 

a Simulation study of the properties of maximum likelihood 
estimators for the parameters of the model is given. 

The correlation structure is derived for other models in 
the GLARMA(p,q) family: the first-order moving average, the 
second-order autoregressive, the first-order mixed autoregres- 
Slve-moving average and a bivariate first-order autoregressive 
. process. These different models, particularly the bivariate 
extension, demonstrate the flexibility of the GLARMA(p,q) 
model. 

In Chapter III the first-order moving average process with 
Exponential marginals of Lawrance and Lewis [Ref. 5] is ex- 
tended to a two parameter model. This is done by utilizing 
the NEAR(1) structure which combines two independent Exponen- 


tial random variables into a random variable with Exponential 
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eochioution. A fairly complete set of characteristics of 
this model are derived. In particular the correlation struc- 
ture, the quantity P(X 417 ¥y) the Laplace transform of 

sums of xX Ss the Laplace transforms of the distribution of 
eeumcs, the (Bartlett) spectrum of counts, and the joint 


Laplace-Stieltjes transforms of X_ and Xn are addressed. 


+1 
These characteristics are compared to those of other proc- 
esses which produce marginally Exponential random variables. 
In Chapter IV the models of Chapter II are used ina 

preliminary data analysis of wind speed data. This repre- 
sents the first effort to apply these models to a large, real 
world data base. A model for simulating wind data is pre- 
sented and parameter estimates for the data are derived using 


the numerical approximation to the maximum likelihood presented 


in Chapter II. 
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IL. MODELS WITH GAMMA MARGINALS 


Pee LNT RODUCTION 

There have been several schemes suggested for the modeling 
of dependent random variables with Gamma marginals. The 
Gamma autoregressive process of order one (GAR(1)) by Gaver 
and Lewis [Ref. 2], the discrete autoregressive process of 
Order one (GDAR(1)) by Jacobs and Lewis [Ref.12], the Gamma 
Beta autoregressive process of order one (GBAR(1)) by Fishman 
[Ref. 13] and Lawrance and Lewis [Ref. 14], and the Gamma auto- 
regressive process of order one (GLAR(1)) by Lewis [Ref. 10]. 
There is also an attempt to use multivariate Gammas obtained 
by the inverse probability integral transform in a time series 
context by Schmeiser [Ref. 15]. 

The GAR(1) model generates an Ce series using the stan- 
dard first-order autoregression equation (first-order stochas- 


tic difference equation) 


- oo OTe CE Pea) 


The innovative factor, Ens has Laplace-Stieltjes transform 
H and the (es random variables have marginal 
Gamma distributions with shape parameter k and scale parameter 


\. The marginal density function of the {X,} random variable 


is 
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fy (eis, A) = £..(x7k,) = x @ Ax (Tree) 


ee mpece >. Ormk: =: 0:, 


The model tends to produce runs of decreasing value when inno- 
vative term has successive realizations of value zero. The 
GAR(1) is in this sense highly degenerate, even though it is 
a true linear process. Ad hoc estimates of model parameters 
are available which produce the exact 9 value if the series 
is long enough [Ref. 2]. However, maximum likelihood esti- 
mates have not been produced. This model is not extendable 
to a moving average process. 

The GDAR(1) produces an il sequence using the first- 


order autoregressive equation with random coefficients. 


Ko eX _| Se lsvaic., (Gite entee) 
where OV n=1,2,...} is an iid sequence of binary random 
waeraples with EMG) - Tes2 (ON 10) = 0, UG eee Re eee 


- 1S an 1id Gamma sequence. 

This sequence produces runs of constant value when suc- 
cessive realizations for Ve produce value l. When Vn, equals 
zero, a new value is selected. Obviously, this model has 
limited value in general applications and is even more degener- 
Gee ean GAR(1l) process. 

The GBAR(1) is the most flexible model in that it contains 


the GAR(1) and GLAR(1) models as special cases. It produces 
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an mes sequence using 


X = BBX cy Giana 4) 
where {Bos n= 1,2,...} is an iid Beta(k-q,q) sequence. ¢ 
was shown by Lawrance and Lewis to be the sum of a Gamma 
variable and the innovation process of the GAR(1) process 
[Ref. 14]. Although flexible in the sense that it contains the 
other models, it can not be extended to a moving average proc- 
ess. In addition, conditional densities and, hence, maximum 
likelihood estimates are not available. This is because the 
innovation random variable for the GAR(1l) process, while it 
can be generated as a random sum of random variables, does not 
have a known distribution function. 

The most valuable and flexible model seems to be the 
GLAR(1) which produces an ie sequence using the stochastic 


difference equation with random coefficients 


career tn Gay Giaaeres), 
where Oe, n= 0,l,...} is a second-order stationary sequence 
of Gamma random variables, {By ee 1 iss tae {Cy T=, lees} 


are iid Beta random variables, and Gr is an iid sequence of 
Gamma random variables. This model is extendable as an auto- 
regressive process of arbitrary order and as a moving average 


process of arbitrary order. These two forms can also be combined 
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to form a mixed model, the so-called Gamma Lewis autoregressive - 
moving average process of orders p and q (GLARMA(p,q)). 

This chapter of the thesis examines some of the special 
cases of the model. One case in particular, the AR(1) form, 
1s reasonably extensively examined. The correlation struc- 
ture is developed. The conditional expectation and density 
are derived. The latter is used as the basis of a numerical 
approximation to the maximum likelihood method of parameter 
estimation. Directional moments and the probability of x +l 
being greater than XO are derived. Ina later Chapter of 
this thesis, this model is used as a basis for analyzing wind 
speed data. 

The special case of the moving average of order one is 
examined in seme detail. The correlation structure is derived 
with some emphasis on exploring the restrictions on the range 
of correlations that are possible. Directional moments and 
an empirical examination of the probability that | 1s 
greater than Xx are examined. 

As a demonstration of the flexibility and extendability 
-of the model the mixed model of order one, the autoregressive 
model of order two, and a bivariate model are introduced and 


their correlation structures derived. 


B. FIRST-ORDER AUTOREGRESSIVE BETA-GAMMA MODEL, GLAR(1) 


ieee Eni TOoduct Lon 
The first-order autoregresSive Beta-Gamma model is a 


special case of the GLARMA(p,q) model when q = 0 and p = l. 
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The autoregressive model generates an {Xx} sequence using 


r BaEn! (I1.B.i.1) 


where ix, n= 0,1l,...} is a second-order stationary sequence 


of random variables with a Gamma(k,1l) marginal distribution; 


{AL, n= 1,2,...} is an iid sequence of Beta(k-q,q) random 
variables; (BL, n= 1,2,...} is an iid sequence of Beta(aq,k-q) 
random variables; {G_, m = 071,.:..) is an iid sequence of 


Gamma (k,l) random variables; HENRI ee and {G} are inde- 
Pemeemt; 09 <q< k. 
Choosing X 


= G, makes the {X,} sequence stationary. 


0 0 
The Gamma random variables are parameterized by the shape 
parameter and the mean, rather than the scale parameter. MThis 


somewhat unusual parameterization has some advantages in 


statistical work since Gamma(k,u) = u Gamma(k,1l) [Ref. 7]. 

A Gamma(k,1) random variable has density 
£_(x:k,l1) = ee es, (SET BD 
Cn Se ene) ’ t etre ie 


This is a special case of the more general density 





(Kk _Kx 
: eel” k-1 U 
£. (sik, u) Fiky * e 
where 1p = 1. Of course, since the scale parameter, shape 
parameter and mean are related by the relation u = <, the den- 


sity can be specified by any two of the three parameters. The 


PAF 





typical parameterization in terms of the scale parameter, A, 
is useful because G(k, ,A) + G(k,,4) = G(k, +k,,A). This re- 
lationship is not true when the parameterization is through 
the shape parameter, k, and the mean, u. 


A Beta(m,n) random variable has density 


=  ftmrn)_ ,m-l 4.,,n-1 
f, (xim,n) = T(m)I(n) x (i x) , (leer so) 
Oka eR ae! re (¢ SOU Sn O oe Gem O58 


For the Beta random variable to be properly defined 
each of the parameters must be positive. Hence, when g = k, 
(II.B.1.1), as defined above, is no longer appropriate since 
each Beta random variable has a parameter that is identically 
zero. In this case when q = k it is understood that the {A_]} 
sequence is considered to be identically zero and the {BL} 
sequence one. Therefore, II.B.1.1. becomes simply x = G, 
and the {x} sequence, like the come sequence, is iid. A 
justification of this generation scheme for a Gamma process 
as defined by II.B.1.1 was provided by Lawrance and Lewis 
(eer. 14, pp.24). 

In this section the correlation structures of the oe 
sequence and that of the oe and {G3 Sequences are addressed. 
Other characteristics of the sequence, such as conditional 
expectation of x given Rea directional moments, and 
P(X, >X) are considered. Of particular note is the deri- 


+] 


vation of the conditional density of X, given X,_1- This 
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leads to the formulation of the likelihood function and a 
computer program to generate maximum likelihood estimates 
of parameter values. The numerical convergence properties 
of the likelihood method are assessed in Section II.E. 
meeecOorrelation Structure 
The serial correlation of the ee series is easily 


determined by a straightforward calculation. We have 


Now X. and G. are independent if j > i and Xs and A. are 
independent if j > i. Using these facts along with the iid 
nature and independence of the {BL} and Ga! sequences yields 


the following expression when expectations are taken. 


2 
B(X,X,,) = B(A,)E(X¢_,) + B(B,)B(G)E(X,_4) 
. (ea. RUN a ad 
e 2 k 
_ k7tk-kg-gtgk _  k*tk-g 
Ke eo 
Therefore, 
AHI 
n’*n-1 Hi “2 
and 


Ze 





CORR(X_,X_ 0<q<k. (Tame 2.1) 


=i k 
pis is consistent with the fact that for q = k, {x} isa 
sequence of 1.i.d. Gamma variates which implies that 
CORR(X,X__,) = 0. This correlation is easly extended by an 


induction argument to yield 


Se eres 
CORR(X_,X_) ce) nh 2m? 0. Cee ar 


aie 


The two sequences {x} and Mens Can be viewed as a 
bivariate pair (XG) of Gamma(k,1l) random variables. 
Therefore, the correlation structure of these sequences may 


be of interest. Proceeding as before 


G 
mn A nn-l-n el el 


Taking expectations as before 


2 
E(X_G_) - EB (AL) E(X,_j) E(G,) + E(BL)E(G_) 
meme) = “24.7.1 q.k(kt+l) _. keg , Kgtg 
nn x4 k «2 k “2 
- + 
mec )oo= -— GrimePe2.3) 
nen Ke? 
Therefore 
pric 
COV (X_,G_) 2 


and 
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— g. 
CORR (X_,G_) eo (rt. Be 2.4) 


When q = k this is 1 since x = Gee 


Pursuing the process one more step we determine 


x hs, = ov el ne i Saal 


Taking expectations as before and using II.B.2.3 and the second- 


order stationarity of the {x} sequence, we get 


E (XG. _}) E(A_)E(X,_jG,_)) = E(BL)E(G_)E(G__)) 
2 3 2 2 2 
= (SSS SEE ok —~ ki tkg-k“g-g +gk 
k Z k 3 
k kK 
3 2 
— k +kg-q 
5 
k 
Therefore, 
COV (X_,G ) = q (k=q) (TB ae) 
n’-n-l ‘> 
and 
BERR Coch = fee Ein i a) 
n’-n-l k k 


imimaee.s through 11.B.2.6 can be used in a simple induction 


argument to yield the general result 
ao) (aay = ese TNs 27 
CORR(X_,G__.) (x) 6 ” ), ™m 0,1, ee ( ) 
When j is greater than 1, X. and G. are independent. Hence, 


CORR(X;,G,) = 0, j > i. 
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Breereonal tional Expectation and Conditional Densi cy 


The conditional expectation of X, given X,_, = y is 
trivially determined from the defining equation and the 


moments of the Beta distribution as 


E(x,|X,_p=y) = Aby+ 2. (II.B.3.1) 


Recognizing sd as the correlation between x and X and 


n-1 
Hetting po be that correlation, II.B.3.1 can be written as 


~LQ 
l 


B(X,|X,_ =y) =) fois = poy + (l-p). (LEG. Soe) 


Thus the regression is linear in y. 
The conditional density of x given ET can also be 
determined. It is easiest to start by deriving the condi- 


erenal distribution of x given Gers = y. 


eo a) 060 CPA y + BaGec x) a= P(BeG. <x -A_y) 


Now [Ref. 14] the product of a Beta(q,k-q) random variable and 
a Gamma(k,1l) random variable is a random variable with density 
Gamma (q, 4) . Hence, if we let dD. be a Gamma (q,#) random 
variable, 


P(X, < x|X eS PD ae CED Bois 2) 


n-l 7” 


imawsecan be written as a convolution if one is careful about 
the upper limit of integration. Since Do is Gamma (q,#) , it 


is non-negative. Hence, PDE aA) TSeZ2ero. ie x-ALY is 


o) 





less than zero. Since AL is a Beta random variable and, 
hence, bounded above by one, x-ALY can not be negative if 

x > y. However, if y > x, then A, must be restricted to lie 
in the range (0,5). Taking this restriction into account 


and writing the RHS of II.B.3.3 as a convolution 


E 
Bee = x) X__, =) = J £ (2) F,(x-yz) dz, (TT Bes) 
and 
Sooo See een gr 
a= - 
— nex <7, 
y y 


where £, (2) is the density of a Beta random variable as in 
Mees. 1.3, and Fy is the distribution function of a Gamma 
ge random variable. 

GEecOuLse, to get the conditional density of xX 
given Xl = y, we must take the derivative of II.B.3.4 
with respect to x. Recognizing that the upper limit may 


be a function of x and applying Leibneitz's Rule where 


maperoprrate yields 


L 
aS x ey ee Cee aye, Aye (Ee Begs) 
n'“n-1 0 
and 
al de xX > Y, 
oF = 
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where f (x) is the conditional density of xX iven 
‘| n J 


X-y! £. (2) is the Beta density, and fy (x-y2) is the Gamma 
(q,¢) density as in II.B.1.2. Writing this result in terms 


of the densities involved we have 


L 





Z r' (k) k-q-l,,__\q-l_k? qg-1_-k(x-yz) 
rex, J T(x-q)T (aq) le 2) Tae e dy, 


Meth the condition on L, as in II.B.3.5. And finaly 


Ee Lx (x) = {1 5} x Fe“ 
n'~n-l LU (ke=oeatl Ce) 
et a -~1k 
See IT" (1-2) PT" (x-yz) Te Yay, (II.B.3.6) 
0 
and 


As a check on the derivation of the conditional density, 
*the conditional density and conditional expectation were calcu- 
lated for values of k and g which produced simple integrands. 


One of these cases was k = 2, q= 1. Then k-q-l = 0 and 


q-l = 0. In these parameter values II.B.3.6 reduces to 
13 feo -= 2 if Y4qaz Lb = 
n= i 0 | x 
7 itn < << ¥ 
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Pecet integration, 





2y 
= ae a 
—— EEK ye 
[ 7 7 y 
Bai) (II.B.3.7) 
n'n-l -2x 
ee tt x < y 
iy Y 


Since the two expressions in II.B.3.7 are non-negative, we 
can insure that this is a density by verifying that its inte- 


feeal 1S one. 








a fy jy (dx = | 2 Saray ote eine 
0 n'°n-l of y vy 
y y 2 = ae 
= i J ax - i f e 2Xqy - = - iyi e 2Xa, 
= ae one ae _— oes 
2y 2y 2y y 
= eg 


as required. We can also take the conditional expectation to 


see if it equals s + = asmuequneed by TlLabe3 . barbus, 


% v4 i a = “ -2x aoe IL 
= = - ——J]dx + xe [(—— - —]dx 
(x) dx J ele = ] I 7 y 


bes) 











J BE GK Ste ee OR xe ?*ax i — i —2x 
: eno f f [ 7 =. Sa ax 
«Oe ee Ge ee 
7 Z 4y Ay 2 4y 
7 ers _ ane 
Z Ay 
a Greg 
eo 
as required. 
Using k = 3 and q = 1 produces a density of 
jaro (e3¥-8 -- 25 if > 
y oye Syn eee 
rt [ox (x) = (dig ys ee Sore oy), 
me n—L 
2x 2 =k 
= - — (1 - ) Wee ee yee 
2 2 
ny 3y 
This density is non-negative and integrates to one. It also 
produces a conditional expectation of Y + 5 as required by 


1010s ieee eee 
results obtained from numerical integrations of II.B. 


~ estimation applications. 


These results are also of use in validating the 


B62 a 


This conditional density can be used to form the joint 


density of Ky rXoree- 


of the data. 
section. 


4, Maximum Likelihood Estimation 


1X, and, hence, the likelihood function 


This subject will be addressed in the following 


Once the conditional density of X_ given X,_, and the 


marginal density of X, are known, it is possible to evaluate 
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the joint density of Ky rKore+ +s X- Since the first-order 
autoregressive process is Markovian, as can be seen directly 
from the defining equation II.B.1.1, the following equation 


is valid: 


RS, Xe cp tw yp &: = 5 
( 3 | 4-1 %}-2 x, ) £(x5|x,_)), Wes See (aaa de) 


Applying II.B.4.1 n-l times to the joint density of XX 


vee Ky produces 


E(x Xpiypr eee X) = £ ix] xi) £1 (x1 xp) «+ £, (xa |x) £5 (xy) 


Cities <4. 2) 


where £ is the joint density of Koreee Kyi ay is the condi- 


tional density of x. given he and f. is the marginal den- 


1! 
sity of the {xX} sequence. 
Viewing the joint density as the likelihood function, 


letting L be the likelihood function, and taking natural 


logarithms of each side of II.B.4.2 produces 


in L = In £, (x) + ) ln Sy ce ng [ee CE 4) 


Recall that in Section II.B.3 the conditional density 


was determined to be 


ay 


i 
fely) = (— A pte he 0 Rea-2 gy) TL yz) ky gz, 
Wiek=¢) [P(q) } 0 


| 1 it moe Se (G5 234 24) 


For a given set of data the likelihood can be viewed as a 


function of the parameters k and gq. Let 


= JL 
G(k,q) = In £5(x)) + a In £,(x,; 


$1 |X) (II.B.4.5) 
Assume for the moment that a procedure has been established 
to evaluate G(k,q) given k,g (0 < q < k) and the data. Consider 
eae problem of constructing a program to determine the values 
Of k and q which maximize G(k,q). An outline of the program 
can be constructed as follows: 

1. Read the initial values for k and q. 

Read the data. 
2. Determine a direction of search. Use the following 


difference equations to approximate derivatives. 


; = GUN Sei ele eo 
3 ~ G(k,gtdg) - G(k,q) | 
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mec i ktik, dd) = G(k,G) . Glk,gtAg) - G(k,q) 
ec Dk, mee and Dd; = 


DK), then YG; 


Dd. 
1 
approximates the gradient of G at the current k,q values. 


ier the ith iteration. If we define VG; = ( 


por i = 1, let the direction of search, d be VG. . FOr 


ae 
1 > 1, define the direction of search, qe to be 


Slyeell <i ec Cigtan ore wale, (Gar SB 458) 
VG, VG 

where 8, _4 is —-——__, cic BatetCwOr erie  Fengenmer scilc 
Lie cial 


present and preceding gradients. Formula II.B.4.8 is the key 
equation in implementing the Fletcher-Reeves Conjugate Gradi- 


ent Algorithm. Once d, has been selected, normalize its 


k 
length to one. 
ee Let the initial step length, Sh, be on and let N= 1. 


@Semouce the trial values of k and q, Tk and Tq, using 


Hie k + She DK. , 
Gm. B «47. 9) 
tem = at SL * Dq, - 


imeeG(Tk,Ta) > G(k,q), let SL = 2 * SL, otherwise go to 4. 
fee = 105°>k = Tk, gq = Ta, go to 2. (No step larger than 


Po * 10 > ~1.0 Mena ieoweon Sl ian 0 a Nea oN) ge ter 3. 
4. If N> 1 ( at least one step produced an increase), 


use a golden section search between Rk, leew ror -steop eN—2 vand 


Shy, 





Tk and Tq for step N to determine the maximum function 
value and the k and q values, kMAX and qMAX, which produced 
it. Here k = kMAX, q = qMAX, go to 2. 

ieee = 1, go to 5. 

9. Since the initial step along the indicated direction 
produced a decrease in the function value, check to see if 
you are at a local maximum. Determine the function value 
at points at 30° intervals on the circumference of circles 
Meee radii of em PO On (O° is parallel to the g axis). 
If the maximum of these test values is greater than the pre- 
sent value, set k and q to the values which produced the 
maximum value, set the present function value to the maxi- 
mum, set 1 = 1 and go to 2. Otherwise terminate. 

The program above assumed that, given gq, k and the 
data, the value of the likelihood could be calculated. One 
difficulty in performing the calculation is that the integrand 
of the conditional density may contain singularities. AS a 
precondition to using an IMSL routine to evaluate the inte- 
gral, these potential singularities must be removed. The 
technique employed requires that the coefficient of the term 
that goesto infinity as one of the limits of integration is 
approached is added and subtracted. The part that is added 
is then integrated separately and added to the part that is 
evaluated by the IMSL routine. 

To procede with this technique we first spiit the 


Mieegral into two parts. Thus, ignoring the part of I1.B.4.4 
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outside the integral, we have 


L 
x (xey2) T+ek¥2a2 4 i emcee (l=) (ia ee ze ere a 
Ey 2 


(EES e420) 


In the first part of the RHS of II.B.4.10 the term a Could 
tend to infinity as z tends to zero if k-q-1 < 0. If we set 
z equal to zero in the remainder of the integrand we get 


(1-2) I+ (x-yz) Pek? 35805. Adding and subtacting this 


= _ a 
f 2* q 1 i-z)4 b ixeyz)4 ve aay 
0 
be geal qui So iy eee (ka etic =dek=e=8 
= f Base = (1 =2) (x-yz) e -x1 “z +X Z dz 
@) 
2 - 7 - aio pee 
= gk-a-1) (y-2)4 1 x-y2)? Loky2_,4 Ly az + x4 1 f see Sag 
0 0 
Thus, 
L/2 : 7 
f gel (429 tb ix-yz) 1 KY2 a, 
0 
Da k—G 
yaw epee. = = = = =a) (5) 
=] gl elena t x-yz) 2 tek¥%_ x4 *]dz+x4 9 
0 
(1teBoe 2 i) 
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Recall that since q < k, k-q > 0 and the second integral on 
the RHS can, in fact, be integrated as shown. Now the inte- 
gral on the RHS can be evaluated by IMSL routine DCADRE and 
the second part of the RHS can be easily computed. 
Applying this technique to the second half of the 
x 


integral, recalling that the case where L = 1 and L = 5 


must be considered separately produces 


Ar 
a9 4d (gig) Fh gey zy Pek Y2 az 

iy 2 

' k-q-1 -l kyz q-1_ky gq-l1 

= f i. (x-yz) 7 ey - (x-y) ey | Gla) maiz (Eland col 2) 
Mg/2 
7 (b/s 
+ (x-y)4 TeX q 

and 


5) : = 
2X q * (lara) 2 1 x-yz) 4 T KY2 az 
ey 2 


x/Y a Be = ees =e 
J (x-yz)? 112k q ae 2) 2 TeKY2_ (2) * q "(1-34 ae * dz 
my ZY. ri 
(=) 4 
(Grr 25.47, 13) 
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Since q > 0, the second terms in the two previous expressions 


eeae properly integrated. 
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Two points should be noted. First, if k-q-l > 0 or 
if q-l > 0, then these steps are not required. But whether 
they are required or not, they are always accurate. Since 
the exact path of the search algorithm is unknown at the 
start, these expressions are used throughout to insure accu- 
rate calculations regardless of the k and q values encoun- 
tered. Second, if x = y, two parts of the integrand simul- 
taneously tend to infinity and this procedure breaks down. 
This does not pose a problem for continuous data since the 
probability of this occurring is zero. However, if discrete 
values are used or if the data is truncated because of limits 
on the accuracy of the measuring device, then the data may 
have to be preprocessed to insure that successive values are 
not equal by adding a small increment to one of the values. 

When the program was written, its accuracy was veri- 
fied by three checkes. The case k = q implies independence 
in the basic model since as q tends to k the probability that 
AL equals zero tends to one and the probability that BL 
equals one tends to one. Hence, in the limit , aaa and G. 
» is an iid sequence. The logarithm of the likelihood function 
for independence was calculated and compared to the program 
results for several values of k and q where k = q. The two 
calculations were equal within machine roundoff and compu- 
tational accuracy. The special cases of k = 2, q = 1 and 
[= 8o, g = 1 discussed in II.B.3 were also computed. The 


logarithm of the likelihood function was computed using each 


a5 





of the conditional densities derived in II.B.3.7 and II.B.3.8. 
When the results of these calculations were compared to the 
program results with the specified k and q values, the re- 
sults were equal within calculation and roundoff accuracy. 

The results of the tests of this program when used 
with simulated data with known parameter values are presented 
mamoeceion IL.E. 

Note that there are natural moment estimators for the 


three parameters, k and q and y in this model. These follow 


from the fact that 





= eel 
9 (1) CORR(X_X, 44) LS ee 
Var (X<.) a 
2 eee al Var (X) an et 
2 a a ee ee ane 
[E(X_) ] [E(X) ] L 
Thus we use for moment estimations 
( oes (II.B.4.14) 
qg = (1-p(1)k Gere) 
nw Gz) - 
~ = “2 (II .B.4.16) 
S 
x 


These moment estimations will be compared to maximum likeli- 


hood estimations in the case where u = 1 in Section I1I.E. 
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5. Directional Moments 

Unlike processes with normal marginals, non-normal 
processes are not completely determined by their correlation 
structure. Directional moments not only demonstrate this 
difference, but also help to differentiate processes with 
Similar correlation structure and identical marginal distri- 
butions. They can also be used to help determine parameter 
values. In addition, they may also be viewed as another 
mayeoL Characterizing the joint distribution of the process. 
With 
+ BnG_, 


io n° n-l 


with all random variables defined as in II.B.1.1, we first 


address E(X x* ). We have 
Mm n—L 
Xo Axn-1 + B ci 
2 = 3 2 
el AOS aren a 


Taking expectations, recalling G; and Ra are independent if 
> 4, JES and acer are independent, and A; and x. are 


independent if i> j. Then 


i 3 Zs 
By(KeX Se) E(A,)E(X>_,) + E(B,)B(G,)E(X2_)) 


k-q k(k+1l) (k+2) , g , 1 . K(k+) 
k = “ ee 


nw 
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1? 5 


2 _ k(k+l1) , (k-q) 
BUX Xn) = ARS ED | cee2) + a 


2 
2 k (k+1) -k° +2 (k- 
Beex= 5) = ao ones eo (II.B.5.1) 
Since 

2 ae. iO: 
Ky = oil - eo, 1 iu Bon! 

2 eo: 2 Do: 
sn 7 Ad a AB Sn Xn-1 7 ee nol 


Taking expectations as before 


2 2 3 2 
E(X°X__,) E(AD)E(X7_)) + 2E(A,)E(B_)E(G_)E(X~_,) 


2 2 
+ E(BY)E(G™)E(X,_)) 


(k-g) (k-gtl) k(k+l) (k+2) , 2(k-q) gq.) k(k+l) 
iS (ican) x? k k ye 


in eter) K(k+l) 5 
(ic ee) eer ec 


After simplification we get 


2 = Seite Ges meee cee) B.S. 
E(X-X_y) — so Cis Pgou) 


Note that these two directional moments are different, 


indicating that the process is not time reversible. 
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PIQtnNer Characterization Of the JOINnewarserroution 
and the time directionality of the process is given by 
> xX). There is a simple analytical solution for 
meee in the GLAR(1) process. Consider, 
ee! 


>X ) = P(A Xx 


n+l n n+1“n aril meee a 


= P(B > [1-A_ 3 /X,) GEL .5 Geel) 


ee ee 


Recall that BL Beta(q,k-q) and ss real is Betatkea, as. 


yi 
Hence, [1-A_ 4] 1s Beta(q,k-q). Since Co and x are inde- 
pendent and have the same marginal distribution and SA aeay and 


B are independent, each side of the inequality in II.B.6.1 


n+l 
has the same distribution and the random variables are inde- 


pendent. Hence 


This is a strong property of the process. While the 
process, as seen by the directional moments, is not time 
reversible, the fact that Pee) = 0.50 means that 
sample paths will have as many "runs up" as "runs down". 
Sample paths are given in Figures II.B.6.1 and I1I.B.6.2. 

An additional point of interest occurs when k = l 


and the process has exponential marginal distributions. 
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Another exponential first-order autoregressive process in 
which P(X 7 X) = 0.50 is the PREAR(1) process. This is 

a special case of the two parameter NEAR(1) exponential 
process of Lawrance and Lewis [Ref. 8] in which the two 
parameters a and 8 are related by 8 = —. The two 
exponential processes are very similar in sample path proper- 
ties. However, the GLAR(1) process has a smoother joint 
Seseribution. In fact, the likelihood for the PREAR(1) 


process is discontinuous, making it difficult to get maximum 


likelihood estimators. 


C. FIRST-ORDER MOVING AVERAGE BETA-GAMMA MODEL, GLMA(1) 
ime Introduction 
Another special case of the GLARMA(p,q) model is the 
first-order moving average model where p = 0 andqe=4l1. This 
arises naturally from the key result that an tX} sequence 
can be formed by the random sum of two independent Gamma ran- 
dom variables. In the first-order autoregressive case of 


Section II.B, the generation scheme was given by equation 


II.B.1.1 and is repeated here 


Mae distribution of the Oke sequence depends on the inde- 
Pemaence and distribution characteristics of ae and G.- 
It should be noted, however, that any two independent random 


variables with the required Gamma distribution could be 


50 








substituted for Xl and Ge without changing the distribution 


and G 


Opa Xo: Mirmpartieculdr if we substicute Go OE x 


aa 


mor Gir we produce the first-order moving average process 


n-l 


which generates the es sequence using 


ry rs) oe La 


where oY, n 1,2,...} is a second-order stationary sequence 


ereecancdom variables with marginal Gamma(k,1) distributions, 


{AL, n= 1,2,...} is an iid sequence of Beta (k-q,q) random 
variables; {BL hee 2 Swans UdmucenqMence Of Beran G co) 
random variables; {G_, = 07 ees wane id ecoguemecre £ 


Gamma (k,1) random variables; {AL}, LBA, and ices are inde- 
pendent; 0 <q< k. The Gamma random variables are para- 
Mmereerized aS in II.B.1 with density as in II.B.1.2. The 
Beta random variables have density as in II.B.1.3. 

In this section we will address the correlation struc- 


ture of the {X,} sequences and that of the {X J, {G} se- 


Guences, theoretical ranges for possible correlations for 


the eS sequences, directional moments, and the (OS ay = | 
2. Correlation Structure 
Using II.c.1.1 to define X, and X__, we have 
[oon-l (A,Sy 7 ene nena a ere one. 
2 
ee eonoread © a pein nate ne lence sane n= 1 nee 


Syl 








Using the iid nature and independence of {ALt, {Bot, and 


ace) and taking expectations 


Bex 4) = E(A,JE(A,_j)E(G)E(G,_,)+E(A,)E(B_j)E(G,)E(G__,) 
Zz 
+ E(A,_,)E(BL)E(G_,)+E(BL)E(B._j)E(G__,)E(G,_») 
ee (aye k-q, (q k~-q, ,q, ,k[k+1] 2 
(Ay* + (SA) A) + (Fa) + 
k- 
= 1+ (£4) (g) (2) Gate. 2 I) 
Therefore, 
_ a ei Al 
cOV(X,,X,4) = (Sb) (Dp 
and 
Z mec twicl 
CORR(X_,X__j) (Ob Fy) ae) Gigeey) 7) 


A simple calculation will show that for lags greater 
eiemmeonme the correlation is zero. So equation II.C.2.2 plus 
the knowledge that greater lags are zero is sufficient to 
specify the correlation structure of the {XJ sequence. 

One might note at this juncture that the correlation 


of the {x1 sequence is constrained to lie in the interval 
i 
- 
laxing it will be discussed later. It is also worthy of 


(0,—). The reason for this constraint and a method for re- 


a2 





note that the {X_} sequence is stationary and has the same 
marginal distribution as the es sequence, if Xp = Gp. 
Ase indicated in IT.B.2 the Ue and cGee sequences 
can be considered to be a bivariate, correlated Gamma(k,1) 
process. As such, the correlation structure of this bivari- 


ate Gamma may be of interest. Consequently, we first develop 


mie correlation of x and ee in the Standard fEashion. 


x = ALS _ BASn 1 

2 
Xx G = AG +B GG 
mon mon nnon-l 


Taking expectations as before, 


= 2 
E(X,G_) = E(A,)E(G2) + E(B))E(G,)E(G,_)) 
- (Kogy (KUstll) reac 
k 2 k 
k 
2 
ee a eee | 
7 2 
k 
Therefore, 
_ kg 


and 


3) 6 





k= 
CORR(X,,G_) = 4 ; (jC 3) 


Now consider the correlation of Ka and G. 1° We have 


Expectations in the standard fashion produce 


2 
E(XG__}) E(A_) E(G_) E(G,_y) oF E(B) E(G__y) 
. k-g 4 q Kik+ll) 
k k “2 
2 
EP Sek 
“2 
Hence, 
age 
COV(X_,G__)) 2 
And finally 
a 
CORR (X_,G,_1) = eee Cie sGrenes) 


A simple calculation convinces one that the eorrelativen for 
lags greater than one is zero. In addition, it is clear 


Peenat CORR(X_ Gi.) =lOe fOr tmeowly2 ats, —HeNGey ete. 2 and 
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1.C.2.4 are sufficient to specify the correlation structure 
of the {xX} and Nene sequences. 

It has been noted before that the range of correla- 
tions for the first-order moving average process generated 
Eyeene Becta-Gamma method of II.c.l1.1 is constrained to lie 
in the interval from zero to one-quarter. There may be other 
random linear combinations of Gamma random variables which 
give a moving average process with Gamma marginals and a 
greater range of positive correlations. Thus we now examine 
a more general hypothetical generation process to prove that 
any random, linear combination of two independent Gamma random 
variables which generates a sequence with the same first two 
moments as those Gamma variables has a correlation that lies 
miethis same interval. In fact, this proof only requires 
that the dependent random variable have the same first two 


marginal moments as the generating Gamma variables. 


PnisOREM : 


If the ee) sequence is generated by 


= ie. Cazeo 
ae ALG % ono ( 
where {XJ} is a second-order stationary, non-negative sequence 
of random variables with the same first two moments as the 

ase} sequence; ae and {BL} are iid sequences of random 


variables; {G} is an iid sequence of Gamma(k,1l) random 
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variables; and Ue {B_}, and {G} are independent, then 


0 < CORR(X,,X_,) < 0.25. 


Proof: Since {a}, (BI, and {G_} are independent and ae 


is non-negative, Ce and {BL} must be non-negative. Hence 


Ea) > 0 and E(B) > QO. 


Taking expectations of II.c.2.5, we have 
E(x) = E(A )E(G,) i E(B) E(G 
ie =L FEA) + E(B) Ge Eas Gite ors oi) 


Peemee, O < E(A) < 1 and O < E(B) < l. 


Computing the serial correlation of ioe yields 


+ B G ) 


a4 SeehaGut BaGee MA, Cae n-l“n-2 


non-l bane 


+A 3B ce > Bags 


5 AL An-1enen-1 TA 2 G n-l-n-n-l n ak eine? 


n n-1°n n-2 


Then 
E(XX)_1) - E(A,) E(A,_1) E(G,) E(G,_) + E(A,) E(B, _,) E(G,)E(G__.) 
Z 
~ E(A,_,) E(B,) E(G)_)) + E(B,) E(B, _j) E(G,_1) B(G,_ 9) 


I 


woe E(A) E(B) (z). 
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Since E(x) = E(G_) = l, 


= 1 
COV(X_,X, = E (A) E(B) (=) 


=i’ 


And since VAR(X_) = VAR(G_), 


CORR(X_,X), Se a Ge) 


at 


Using II.C.2.6 and its consequences, we have 


Pee 2 275) 


oe 


Sere a(n =" (A) je 


Se, in general, if Pe is second-order stationary with the 
same first two moments as uGat ERewoerelalmcCorcne lation e£ 
3 is bounded below by zero and above by one-fourth. Q.E.D. 
This constraint on the correlation appears to be 
restrictive since, in the classical case, when two normally 
distributed random variables are added to produce a normal se- 
quence, the range of correlations is (= 545) . The two situa- 
Seems, however, are not comparable. It is clear upon reflection 
that the constraints imposed on the Sa sequence in the pre- 
vious theorem are more severe than those imposed upon the 
classical normal case. In the above theorem we required that 
both the mean and variance of the CO sequence equal that of 


the innovative sequence. However, in the classical normal 


case (where zero mean normals are used as innovative factors) 


Di] 





only the mean of the generated sequence is equal to the mean 
of the innovative sequence. The variances are not equal. 

We now examine the case where the generated sequence is re- 
quired to have the same mean as the innovative sequence, but 
is free to have a different variance. (This 1s the case 
with the usual constant coefficient, linear additive MA(l) 


scheme.) 


THEOREM 
If the non-negative 1X sequence is generated by 
II.C.2.5 and all variables are defined as for that equation 
except that {X,} is only constrained to have the same first 
moment as iG}, ene Oe< CORR(X_,X,_4) eee 
Proof: Taking expectations of II.C.2.5 with the new circum- 


stances produces 


E(X) = £=[E(A) + E(B) JE(G) 


ame@ee < E{A) < 1, 0 < E(B) < l by following reasoning identi- 
cal to that above. Calculation to determine the serial corre- 


Macion Can initially proceed as usual, 


3) 








COV(X,,X,_,) = E(A) E(B) (=): 


To this point all calculations and reasoning have been the 
Same as that which produced II.C.2.7. However, since {x} 
is not constrained to have the same second moment as {G_} 


the most explicit result obtainable is 


E(A) E(B) (z) 


owe VARS 


where VAR(X) is, of course, a function of VAR(A), VAR(B), 
and VAR(G). Since it is obvious that the smaller the value 
for VAR(X), the greater the serial correlation for {xX}, 

let us reduce VAR(X) to its smallest values. Since the dis- 
mimrpution of iG} has been specified, its variance is fixed. 
Let P(A, = a) = 1 and P(BL=b) = 1. Then trivially E(A) = a, 
E(B) = b and VAR(X) = (eee be neo Under these conditions 


k 
II.c.2.8 becomes 


If we further specify that a = b, then CORR (X_,X__y) achieves 
mecemaximum Of 1/2. Q.E.D. 

The situation developed above is comparable to the 
classical situation where the innovative factors have distri- 
beeen, nO pao Except for the degenerate case where one 
coefficient is zero, the sequence generated by a linear com- 


bination of innovative factors may have the same first moment 


a8 


Cenc. 





as the innovative factors, but will have a different second 
moment. So, under comparable conditions the random, linear 
combination of Gamma(k,1l) random variables can produce 
positive correlations equal in magnitude to the positive 
correlation produced by a classical normal process. The 
distribution of the ee under these circumstances is unknown. 

If the distribution of the gi is constrained to 
be Gamma(k,1l) and the Beta-Gamma generation scheme is used 
to generate the {X}, then the maximum correlation that can 
be achieved is one-fourth. 

3. Directional Moments 

As mentioned in II.B.5 the directional moments of a 

non-normal process are not necessarily equal and may provide 


valuable information about a time series. First, consider 


B(x*X eee GON C1. 1 
n n-l 
eek 22 DP 

Xx = AEn + 2A BS Snel ~ BASn-l 
and 

ae = eal n-l i n-l-“n-2 
Therefore, 

2 - 2 Wy Z aS 3 ce 
er 1 7 AD An-12n Gy ae eS ae ee n-ln n-2 
2A BB +B°B. .G* _G 
oe ae n- 1° n G.-1° n- el er richge ho en leas eo 
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Peeing expectations as in II.C.2 yields 


(k-g) *(k-g+l) , 2(k-g)* (k+l) | (k-g) q(t) (k+2) 
3 cS 4 


E(X°X 
‘ i k 


ajay 


. beg) (k-gtl)g , 2k-g)7a* | gM (qth) 
kK? “> K? 


Upon simplification this produces 


2 
B(X°X = Oats p, inane (x=) q (2g (kt1) +kt2, 
n ntl 3 4 k 
k k 
2 
if 4.[2k-qt1] (rene eT) 
k 
In an analogous fashion 
2 = 2 2 2 Z 
to n- 1 . rene lon n= 1 ‘i rT ne ig n- 1°nSn-1° n-2 i eae enn—2 


ia 3 2 2 2 
eo cee * Se eee Cee ene? 7 Oe ned ne lon}? 


Taking expectations we have 


etx x2.) = deeg*Gerg+l) , 2 (k= act ‘aq , (kg) qigtl) 
ie ti L K “> Ke 


2 2 
_ (seg) (e-gt2)g(c+2) , 2Uk-g)q? (k+l) 5 g@(gth) 
KA = 


which simplifies to 


Gull 





2 
2 = Usesegr ie Xe fla) Suigge (2 ioc) ) 
en -1? 3 : 3 


: 2 
nt Ahm gil) @erec3 2) 
k 


a Empirical P(X 417 Xp) 
No analytical procedure was found to determine 


P(X > XX). Hence, a simple computer program was constructed 


n+l 
to evaluate this condition for a series of sixty-eight thou- 
Sand pairs of numbers generated by the Beta-Gamma scheme for 
each of ten random number seeds. The answer obtained was 
considered to be accurate within 0.001. The comparisons were 
run for each of seventy-nine values of q, from 0.05 to 3.25 

in increments of 0.05. All of the results of the comparisons 
fell in the range 0.499 to 0.501. Fourteen of the values were 
different from 0.500. No pattern was apparent in the devia- 
tions from 0.500 and these deviations were considered to be 
Bancom fluctuations within the given margin for error. It 


thus seems clear that ie (Og eee) for this process is like 


the GLAR(1) process but no proof has been found. 


D. OTHER CASES OF THE GLARMA(p,q) MODEL 
i Introduction 
A primary advantage of the GLARMA(p,q) model is the 
ease with which it can be adopted to cover a variety of 
Special cases. Two special cases, the first-order autoregres- 


Sive GLAR(1) and the first-order moving average GLMA(1l), were 
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Severed in I1.B and II.C. The intention here is to briefly 
present three additional cases of the general model and derive 
the correlation structure of each case. The special cases 
considered are the first-order mixed model, GARMA(1,1), 
the second order autoregressive model GAR(2), and a bivariate, 
first-order, autoregressive model BGAR(l1). The purpose 
in presenting these cases is to demonstrate the flexibility 
of GLARMA(p,q) and not to present a complete, detailed dis- 
cussion of each model. Further extensions of the special 
cases of the GLARMA(p,q) model from the examples given are 
obvious. Details will not be given. 
2. GLARMA(1,1) 
Consider the following scheme for generating an ce 


sequence of random variables. 


= ES Egy Oh oan 
x Be An-1 + @nSn’ ( D ) 
= : Tt Lape 2 3.2 
An DAn-1 7 Fanon ( ) 
- where {Xs n= 1,2,...} is a second-order stationary sequence 
of Gamma(k,1) random variables; {AL Tie Ocvilee eee) cesta 


second-order stationary sequence of Gamma (k,l) random 
variables; iB, n= 1,2,..-.} is an iid sequence of Beta(k-q,q) 
random variables; {C, n= 1,2,...} is an iid sequence of 

Beta (q,k-q) random variables; {Dos Newel ee tS an ddece— 


guence of Beta(k-r,r) random variables; Sey Te ele ee eS 
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an iid sequence of Beta(r,k-r) random variables; (Bt, {Cc}, 
{D.}, | ae and Kee are mutually independent; 0 < q < k; 
O- r < k. The Beta(m,n) density is given in II.B.1.3; the 
Sammat(k,]) density in I1.B.1.2. 

Before the serial correlation of the ee sequence 
can be determined, the serial correlation structure of the {A} 
sequence and the correlation structure between the {A,} and 
mG! sequences must be derived. Proceding first with the 


serial correlation of the {AL} sequence, from II.D.2.2 we 


have 
A = DA + FG 
n nn-l heen 
So 
A_A See +FGA 
nn-l n-1l nnn-l 


Werng the iid nature of {Di}, ee) {Gti noticing that when 

mej, OD. and A., F. and A., and G. and A. are independent; 
- 5 ~ J - J 

“recalling the independence of (Fj, Gee; and taking expec- 


tations yields 


eg We ictal) r 
Nee ee 
EAA ) = eo —C——" (Pas. 26.3) 
nn-l “2 
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Therefore, 


_ ker 
COV(A,,A,_1) = == 
k 
and 
CORR(A ,A.) = Kz (ie Dace) 
nN n= k e e e e 


mm fact Ee is just a GLAR(1) process, so the result is not 
mmeeerstng. Using 11.D.2.2, II.D.2.3, and an induction argu- 
ment leads to the general m-step correlation formula 


= k=-r,m 
CORR(A_,AL_ ) = Sa 4 n a m ea QO. (Li. Dive eo 


The correlation structure between {A} and Gee can 
be derived in a similar fashion. However, it is more direct 
to note that since the eit, sequence is the same as the GLAR(1) 
process of Section II.B, the tAt, WG correlation structure 
will be the same as that derived for the (es see, sequences 
mm £L.B.2. Hence, 
16 


_ k=-r,m 
oreo Gae)) = sp! t=p-) 7 on 


Dae 
= GE Be ale) 


iv 
J 
|v 
oO 


meecourse, 1f ) > 1, then CORR(A,,G.) = 0 
Now the serial correlation for the {x} sequence can 


memcrounad. From II.D.2.1 
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XX — (Baa +C 


n Rot ence 


eal we Sono 


= B_B l= eS 


n n-1*n-14n-2 n pel mai ed 


nel mene on 


Using the stationarity of the ca sequence, the iid nature 
and independence of {B i, {Ch, GL the fact that A, is 
independent of G. when j > i, and the fact that {AL} is inde- 
pendent of {Bhi by construction and taking expectations, we 


have 


E(XX,-1) = E(B) E(B (A 5) tE(BL) E(C,_ ) EA ) 


eee eae ee ail 


+ E(B__j)E(C,)E(A,_,) E(G,)+E(C_)E(C_j) E(G,)E(G,_y), 


2 
[a(B_) ]“E(A,_jA,_7) +E (BL) B(C,,_,)E(A,_1¢,_9) 


Z 2 
+ E(B\_,)E(C_)E(A,_») E(G,) +[E(C,) 1“ (E(G,) 1°. 


(Taleb) 


Poem [1..D.2.1 


= [E(B,)E(A,_1) +E(C,) E(G_) ] x 


[E(B _,)E(A,_,)+E(¢__,)E(G _,)1, 
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2 
E(X,)E(X,_)) = [E(B,)1°E(A,_,)E(A,_,)+2E(B)E(C_)E(A_) E(G_) 


Z 2 
+ [E(C_) 1° (E(G) 1”. Celebs 2.2) 
Weang 11.0.2.7 and I1.D.2.8 to compute COV(X, ,X,_4) Vesericis 


A 2 
COV(X,,X,_,) = [E(B)] LE (AL An 2) ~E(AL_}) E(A,_») | 


+ [E(BL)E(C_)][E(A )-E(A,)E(G_) J 


a Shai 


Since UNE, {Gi}, and (2) are all Gamma(k,1), VAR(X_) = VAR (A_) 
= VAR(G_) . Hence, 

_ 2 
CORR(X,,X,_;) = (E(B) ]°CORR(A,,A__,)+[E(B_) E(C,) ]CORR(A_,G_) 


From II.D.2.4 and II.D.2.6 we know this equals 


k-xr 


Bee (eee K=q) (q) (= 
CORR(X_,X_ (Sy S + Fe © (Ones oO) 


1) 


Dicincguitabe2z.. and LE -Db.2.6 in an induction argument 
yields the general m-step correlation of 
Rees Ih (k=g 


2 PA SAG eS r 
SORR(X xX _,) = (4) A) A+) @). 


Recognizing the expression in brackets as CORR(X_,X__j) and 


letting 9 equal this correlation we have 
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_ Kat em | 
CORR (X_,X__ ) = sia CMe ee. ee (LI, AD P20 9) 


Figure II.D.2.1 shows the possible combinations of 
one- and two-step correlations for the GLARMA(1,1) model. 
This concludes the development of the correlation structure 
of the GLARMA(1,1) model. 

3. GLAR(2) 

Jacobs and Lewis [Ref. 12] first developed the follow- 
ing mixture scheme for generating a p-order autoregressive 
processes. We now adopt that scheme for generating a second- 
order autoregressive sequence of random variables. This is 
the special case of GLARMA(p,q) with p = 2 and gq = 0. As such 


it closely resembles the GLAR(1) process. Let 


oes ke + CG, (tapes a) 


where {X, n= -1,0,1,...} is assumed to be a second-order 
Stationary sequence of Gamma(k,1l) random variables; {Boe 
me 1,2,..-.} is an iid sequence of Beta(k-gq,q) random varia- 
‘bles; tC} is an iid sequence of Beta(gq,k-q) random variables; 
me! is an iid sequence of Gamma(k,1l) random variables; {BLi, 
{Clty ica are independent; also a shee aNSL (Sy enka ay P(t. =1) = 
I-P(T, =2) =a. The Gamma(k,1) and Beta(m,n) densities are 
meme in Ti.B.1.2 and Ii.B.1.3, respectively. 

This generation scheme works even though X-] and 
X are dependent random variables. The mixture of | and 


n-2 
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X-2 produced by II.D.3.1 is Gamma distributed and indepen- 
dent of G.- Hence, II.D.3.1 is simply another example of the 
random sum of two independent Gamma random variables producing 
another Gamma random variable. 

Two special cases of II.D.3.1 are as follows. When 
a = 1, the GLAR(2) process reduces to the GLAR(1). When g = k, 
the me sequence is iid. 

The serial correlation of the {XJ sequence can be 
calculated in the usual fashion, assuming stationarity of 


the process we have 


and 


non-l | non-T n-l * a eae 


Using the independence of er and {Git, the fact that X, 


is independent of oes Be and B. when ] > 1 and taking expec- 


> tations we have 


Z = 
Be) 0 = CO E(BL)E(X 4) + (l-o) E(B) E(X,_1X,_5) TE (CEG) E(X,_y) 
= keg. Koike | a ed q 
= oO ren = a ge | JE rey - )E(X,_1%,-9) tree 
Since E(x) = E(G_) iby eacsuip elton Usingsclie Sccond- order 
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stationarity of {X)} 


eee kg) _ a(k-q) (k+1) 
meex 5) (1 = i +a 


meen Simplification 


ak *+ak-akgq-agtka 


en n-1) - k (q+tak-aq) 
Hence, 
1 a (k-q) 
COV(X, ,X,_4) ( lota tka)! 
and 
= ka) 
CORR(X,,X__j) = ata (k-q) 


Ifa= 1, this equals l - re! 
fee tt is clearly non-negative. 


The lag two correlation can be 


Similar fashion. We have 


and 


va 


K? fie —mce chen Lite Ss 2enO. 


(LO sci) 


Since 


calculated ina 





Pec. X 


= k- k- ® 
>) = Mee) E (xen x _ 5) + (l-a) ($4) E(x" _.) ies 


i 


Using the second-order stationarity of the ie sequence we 


can write 


E(X)E(X, 9) = a (SEN E(X, 4) E(X,_2) + (1-0) (AES) (B(x, 9) 1° 
+ Zte(x)1*, 
So that 
COV(X,,X,_5) = a (2%) (E(X,_)X,_5)-E(X,_,) E(X,_))] 
+ (1-0) (EES) (E(x? 4) -( B(x, 5) 371 
+ 2 - Diz(x))* 
= a (SSF) [E(X,_ 1X5) -E(X,_)) E(X,_5) | 
+ fir) (hod) tp (x? 9) - (e(K,_5)}71 
= «(By covix,_,,x,_,) + A varix_,). 
Hence, 
CORR(X,,X,_5) = a (S22) CORR (x, 1'*n-2) + Afro) (erg) 


Ye 





_ k- +ak=20 
CORR(X,,X,_5) = (ey) ea) (Galen 3) 


Miemcooliiten DPrecedure for [1.D.3.3, if followed for 
en! will produce the general recursion equation (Yule-Walker) 
that can be used along with II.D.3.2 and II.D.3.3 to compute 


the m-step correlation. The formula thus produced is 


ae are 
Cree see,) = =~) [aCORR(X, ,X_13_,,) 


+ 


(1-a) CORR(X, ,X Nels (De Eee oeceay 


i 2a 


Gi il ae 


As mentioned in previous sections the (X_,G,) pair 
can be considered to be a correlated, bivariate pair of 
Gamma(k,1l) random variables, Therefore, we proceed to 


derive the correlation structure between these two sequences. 


meem L1.D.3.1 we have 


x, = Bosn-T + C G, 
n 
and 
X_G — oe e G+ec era 
n n-T 91 na 


Recalling the independence of tC and VG and xX. and < 


when j > i, taking expectations yields 


13 





- (k7g Spas K+ 1) 
2 
= Sees 
me 
Hence, 
2 
COV (X_,G_) 5 
k 
and 
So 
CORR(X_,G_) k° Geary DmeS Gae Ss) 
Continuing this process we have 
x, 7 Oe eine " CAS n 
and 
AS n-1 7 See reall v Sana 
Thus 
= k-q _g) (Xz q 
Hence, 
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k- 
COV(X,,G._1) = a(=4)cov(x,_1,G__)) 
And 
= k-g 
CORR(X,,G__,) = a(= 4). 


One further step in this process using an arbitrary 
m-step lag produces the general recursion formula 


k- 
CORR(X,,G__.) = (= 4) [acoRR(x_,c 


n+l-m 


Brguremll.DeS.L displays a plot of the possible com- 
Banations of CORR(X_,X__j) and CORR(X_,X,_5) - Note that when 


a = 1, the GLAR(2) process reduces to the GLAR(1) and 


CORR(X_,X,_>) (AEA) * which defines the lower boundary of the 
_ k-g 

Plot. When a = 0, CORR(X_,X__,) = 0 and CORR(X_,X__5) = fone 

mien Goeseencm zero to one. In the imiteGlOn Tor eile aco, 


-when a does not assume an extreme value, CORR(X_,X__4) does 
not reach a value of one. This is demonstrated by the 


following calculation. From II.D.3.3 


e k-g, ,qtak-2aq, | 
CORR(X_,X__>) (see) lavakeq)- 


iemenis Correlation is to equal one, then 


Us 
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+tak=2a 7 k 
Sahay cq em foc Geiypes <3) 


which quickly reduces to 
q(-2ak-q+2kq) = 0 
Since 0 < q, this requires 
-2ak-qt2aq = 0 
Therefore, 


aL so. 

2 (q-k) 
Since we know that 0 < q < k, a must be less than zero. 
However, a 1S a probability and, hence, is non-negative. 
Therefore, the original requirement in II.D.3.8 cannot be 


Satisfied. Hence, CORR(X_,X cannot equal one. 


n=2’ 

4. BGAR(1), Bivariate Model 

To this point the only examples of bivariate Gamma 

processes presented were those in which the innovation se- 
quence was one part of the bivariate process with the gener- 
ated sequence the other part. The simple random, linear 
structure of the GLAR(1) process makes it easily extendable 
to a variety of bivariate models. We address only the 


Simplest. Consider the following pair of random variables, 


both of which are formed from the same innovation process, (Gi: 


17 





7 aia = Sea teen OC (Gigs. 4.1) 


aa Oe YEG Grr .D 6422) 
where Lee n= 0,1,...} is a second-order stationary sequence 
of Gamma(k,1) random variables; Oe ee = 0 eS acc 


second-order stationary sequence of Gamma(k,l) random varia- 
bles; {B_., n= 1,2,...} is an iid sequence of Beta (k-q,q) 
random variables; tC, n= 1,2,...} is an iid sequence of 
Beta(q,k-q) random variables; {Diy Meee 2 eee Sen a a 
sequence of Beta(k-r,r) random variables; Ue 0 ee oe 
is an iid sequence of Beta(r,k-r) random variable; {Gy 
n= 1,2,...} is an iid sequence of Gamma(k,1l) random 
Variables; {Bl, {Ct, Og elds and Gay are independent; 
Meee < kK; O <q < k. 
This is a special case of a general Situation. In 
a more general case the {x} and ee sequences could have 
separate, correlated innovation sequences instead of sharing 
a single tae sequence. In addition, the {BL} and 1D. 3 
“ sequences and the tc} and Coe sequences could be correlated. 
Before examining the correlation between Oe! and 
{¥}, it will be necessary to address the correlation of each 
of these sequences with the Gas sequence. This is most 
easily handled by recognizing that the relationship of the 
(X,} and ee sequences to the {G_] sequence is exactly the 


Same as the {AL sequence to the {Gf sequence in II.D.2.2. 


ve 











Hence, the correlation structure will be the same. MThere- 


fo alee let = - 
re, we le Py CORR(X_,X__4) and Py CORR(Y_, 


ee 
n=41,2,..., these correlation structures can be written 


without further analysis as 


Ser eeGe) ep = 1 - CORR(X,X__,) = 1-o,? (II.D. 
CORR(X ,G._) = (H(EAY = (1-o,)oy: (II... 
CORR(X,,G._) = (2) (S)™ = (1-p,)o@, n> m> 0. (IIL. 
and 

Peeiye y= = = 1- CORR(Y,,Y,_,) = 1- oy? (IID 
CORR(Y,,G._3) = (F(A) = (1 vy) oy, cra 
CORR(Y.,G 1) = (E) (ALE)™ = (= 90 yeh Sd ee Geers 


.The assumption is that the bivariate process is stationary, 
although it should be noted that starting the univariate 
processes in a stationary mode does not make the bivariate 
process stationary. The initial pair {Xq1Xq} must have the 
bivariate Gamma distribution associated with the stationary 


Markovian process. 


Now we address the cross correlation between ee and 


{Y. Start with II.D.4.1. We have 
719 


A Sh) 


4) 


i) 


.6) 


A) 


ro) 





and 


> Tae 4 = B mp eeGey 


a eal eal henen— 1: 


Using the independence of {C} and Ree and that of Y, and 


G; and C. when j) > 1 and taking expectations 


= B(B)E(X,_1%,-1) + E(C,)E(G_)E(Y__)). 


so that 


Gil. Det) 


and 


Taking expectations as above 


E(X,¥,) = E(D,)E(X,¥,_1) + E(FL)E(G,X)). 
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2 
Using II.D.4.3 we deduce E(G, x) = a4 so that 
| k 


2 
Ger) = (AS) B(x ¥ 4) + QS. (ia,.D. 420) 


Invoking the second-order stationarity of the besa ee | se- 
Guences, using II.D.4.9 and II.D.4.10 and letting w = E(X,Y,) 


and Zz = E(X.Y¥;_,), we have the two equations 


2 
- k-r Se ea 
Zz = (3) w a 2. (me D4 ale) 


MeanG LL.D.4.12 to substitute for z in II.E.4.1li yields 


= 
| 


Z 
amg = 1g Re 
(FS) Sw + G1 + ae ) 


2 2 
ee eee eee kccr), + (koe, tk tae 
We K = 


When EB (XY) Poemsubstituted For w after simplification, this 


produces 
E(X Y = k“q-kart+k rer (II .D.4.13) 
nin = k(kq+kr-qr) ° ees 
Hence 


rq 


COV (X.Y) k (kqtkr-qr) 


Ou 





7 





and 


= See 
CORR(X_,Y_) Reckp ase Meme = G < ky 0 < © < ke 
Clip. 4a) 
_ (l=p ) (1-py} 
(l-o,)+(1-9,,) Py 
(Y=02) (1=0~) 
x = grag ails) 


(i=py) + (=p) By * 


This latter expression follows since for a given k the corre- 
lation structure 1S parameterized by r and gq or equivalently 


by the serial correlations p, and Oy Suveteadt ih. 4. 3 ).ane 


Xx 
mr.D.4.6. 


Viemeanenow SubStitute [L.E:.4.13 into D.E.4.12 to 


solve for E(x ,Y._))- 


-k koa ; 


kq+kr-qr) 


24 0,Gup4 


2 
) psec es a 
none-l k k ( K 


= Ko qtk r+kqr-k“qr-g"r 


k* (kqtkr-ar) 


Hence, 


= kgr-q*r 


COU, Lory) = 
no nai k* (kqt+kr-qr) 


and 


aoe) 
CORR(X_,Y) ) )(——) = CORR(X_Y_) Oy (Ti sped lo) 


= 





Continuing in this vein produces the general formula 


_ 1g k-q,m 
CORR(X,,Y¥,_,) = rastreae oe GU 58 Ly) 


m 
CORR (X__,Y_) Py, on ae On 


Solving for correlations where the eee lags the 
et sequence is similar to the above process, but somewhat 


elebreviated. «Starting with II.D.4.2, 


we have 


Taking expectations as before gives 


may X ) 


nén-1? eee ee ERQGe) G(X 


n-l 


= §(D.)D(¥_)X,_ 


Using the second-order stationarity of {x} and {fy}, 


E(Y X lers isneown £rem [Il.E<4.13, se 
n-l n-l 
Be) | = (Sexy (kogckartk“rtrg, + = 
non-1l 7 ik k {kqt+kr-gr] k 


n?q+k ?r+kgr-k“gr-ar* 


k* (kqtkr-qr) 
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Hence, 


Z 


= kaqr-qr 
cov(Y.,¥__j) = 
k” (kqtkr-qr) 
and 
2 CE See 
CORR(Y,X__}) (Kqtkr-qr) ' ir ) CORR(X,Y_) py. 


Further computations of this nature produce the 


general formula 


2 gr k-r,m 
SoetY Zam! (Kgtkr-qr) | k 


CORR(X,Y,)o, nn >m 


| Vv 


CD eae 10S (Chante Nope ais 


On ieee ke 


To examine the correlation further, note that if 





Py = Py = Os tnenmre..De4.15 yields 
ee 
CORR(X,,Y_) = | Sierats 
fms if p = 0 (i.e., the {x} and me! pE@eesses are iad 


sequences), this correlation is one. For Ce and an to 
be iid, we must have x = Go = Yn? Therefore, this limiting 
correlation does make sense and suggests that perhaps a bi- 


variate Gamma should be used as the innovation process. This 
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would allow for separate control of the serial correlation 
(auto-correlation) and cross-correlation of the Oa, Ehavet ne? 
sequences. 


~~? 


eee tee, XW = X,Y ~ Y~-3)+ the.effect of 
the innovation sequence is slight and the cross-correlation 
between ee] and ee gees weomzZe on La amore complicated 
model than we have addressed here the cross-correlation may 
be controlled by imposing a correlation on the {BL} and cee 
sequences. 
Cross-coupled processes, as discussed in Gaver and 

Lewis [Ref. 2], are possible. These can be used to create 
negative serial correlations in the yo and i) processes. 
E. NUMERICAL CONVERGENCE OF THE MAXIMUM LIKELIHOOD COMPUTER 

PROGRAM AND SIMULATION STUDY OF PROPERTIES OF ESTIMATORS 

Mhe program described in Section I1I.B.4 for computing the 
conditional density function in the GLAR1 process was used 
in two fashions. First, it was tested by using computer 
generated data from a GLAR(1) process with known parameter 
values k, q, and uw. Simultaneously a Simulation of properties 
Sof melee 's K and q for k and q was conducted. Second, it was 
used to estimate the parameters in the GLAR(1) model for real 
wind speed data. Only the first use is covered in this section. 
The second use is addressed in Chapter IV, Preliminary Data 
Analysis. 

Four aspects of the program were addressed in the use of 
the program with simulated data. These were: sensitivity 
of the maximum seeking method to start point, the size of the 


*maximum likelihood estimate 
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standard deviation of the m.l.e and moment estimates of k and 
q produced, and the degree of bias, if any, in the estimates. 
ime addition, normality of the distributions of the estimates 
was investigated using normal plots and Shapiro-Wilks tests. 
Because Of the large computation time involved in obtain- 
mG a m.l.e. the simulation study was small. Three types of 
data generated from GLAR(1) processes were used to exercise the 
program. Each type of data consisted of ten independent sets 
(replications) of 1000data points each. The k and gq values and 
the correlation were varied from one type of data to another. 
Thus the first type of data was generated with ak value of 
4.0 and aq value of 1.0. These parameter values produce a 
correlation of 0.75 (see equation II.B.2.1). The second type 
of data varied the correlation, but retained the same k value. 
Peer or 4.0 and aq of 3.0 produce a correlation of 0.25. These 
values were used to produce data sets of the second type. 
The third type of data returned to a high correlation, but 
used a small k value. The parameter values used to generate 
memo Gata type were a k of 0.75 and a q of 0.1875. These 
mvalwes also produce a correlation of 0.75. Table II.E.1l 


summarizes these cases. In all cases, yu = l. 


DAS hee ce 


CASE k q p 
J 4.0 8, era) 
a 4.0 30 O25 
8 O25 Oe es fe Ores 





*maximum likelihood estimates 
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The Gamma variates were generated by the program LLRANDOMII 
(Ref. 16] and all runs were performed on the NPS IBM/3033 
computer. 

To test the sensitivity of the maximum likelihood computer 
program to the starting point of the search for a maximun, 
each set of data was used in two runs of the program. The 
first run used the actual parameter values k and g as a start 
point. The second run used the moment estimates k and g Of 
the parameter values (see equations II.B.4.15 and II.B.4.16) as 
a start point. The resulting ees parameter estimates k and 
3 were recorded. 

The first case had k = 4.0 and gq = 1.0. The results of 
the computer runs are presented in Table II.E.2 for data of 
the first type. The last column presents the two-dimensional 
distance in the (k,q) plane between the estimates produced by 
the two different start points for each data set. Although 
the values do not differ widely, the relatively large differ- 
ences in some cases indicate that the likelihood function is 
relatively flat near the maximum. 

However, there is another factor which may be contributing 
to differences in final parameter estimates. The calculation 
of the likelihood function for 1000 data points requires the 
numerical evaluation of 999 integrals (see equations II.B.4.2 


pmo ti.8.4.3). The IMSL routine DCADRE was used for this 


evaluation. Two of the parameters in the call to DCADRE are 


* 8 
maximum likelihood estimates. 
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PASeE. 22.8.2 


Results of Search for Maximum Likelihood EStimates 
in a GLAR(1) Model; k = 4.0, gq = 1.0 


Value of 
Data Starting Value Run Time’ Number of Ending Value Maximum Difference 
Set k Se (in Min.) Iterations k a Likelihood (107%) 
4.000 1.000 20 6 3.668 0.882 -214.479 AG 
2.884 0.538 129 i, 3.624 0.886 -214.488 
4.000 1.000 19 2 4.004 0.979 -204.819 4g 
A255 Te Oak 46 5 4.014 0.983 -204.820 
4.000 1.000 104 14 3.451 0.869 -260.332 11 
3.360 02305 IL 14 3.440 0.867 -260.333 
4.000 1.000 16 D 3.977 1,05). 3-206 416 54 
Sn oe 0.947 al is 3.955 1.041 =-206.418 
4.000 1.000 52 5) 4.246 1.232 -304.365 : 
Ae e233 45 6 4. 254llez3>, -—3042366 
4.000 1.000 47 3 4.061 0.998 -211.074 bs 
4.647 Pols 27 4 4.122 1.028 -211.060 
6 4.000 1.000 61 4 3.593 0.891 -229.087 a 
3h, Skehi Oe7o2 iS 9 3.565 0.883 -229.092 
4.000 1.000 oF 5 4.041 0.973 -241.637 a 
Ae IL 1.069 82 7 4.051 0.974 -241.636 
8 4.000 1.000 oF 6 4.024 0.999 -240.432 Bi 
8 CL Svat 0.856 24 8 4.106 1.033 -240.424 
4.000 1.000 SL 3 4.109 1.054 -255.245 a5 
4.398 dee 2 38 5 4 eae eOsl —2e54229 
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the relative and absolute errors allowed for this calculation. 
Practical considerations of computer run time dictated rela- 


: for each of these parameters. 


Mevely large values of 10. 
This error when applied 999 times in the calculation of the 
likelihood function may have lead to the differences in m.l.e. 
parameter estimates. As a test of this hypothesis the data 
set which produced the largest distance between the pairs of 
estimates (data set 8) was rerun with DCADRE error parameters 
set at iL) eee The run which used the actual parameter values 
as a start point ran in 171 minutes and produced estimates of: 
k = 4.102, q = 1.031. The run which used the moment estimates 
as a start point was terminated after 404 minutes CPU time. At 
that point it had estimates of: k = 4.088, q = 1,025. ' The 
G@mecance of 15 x Lo represents a significant reduction in 
the previous distance of 88 x omen It seems from this exam- 
ple that the program can be made less sensitive to the starting 
point by increasing the accuracy with which DCADRE computes 
Pees ineegrals in the likelihood function. Of course, a con- 
Siderable increase in computational cost is incurred. This 
cost is not practical in these simulations or necessary since 
only a rough idea of the properties of the estimates was sought. 
The results of the runs for data type one are also presented 
maeetgure IT.E.1. First, each pair of estimates, (k,q) and (k,q) 
is plotted in the k,q plane. Then each point is projected 


along each axis to more conveniently reflect the marginal 


Variation in the estimates for each parameter. 
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When k = 4.0 and q = 1.0, the method of moments produces 
estimates for k with a sample mean of 3.94 and a sample 
meamcard deviation of 0.55. The statistics for the estimates 
of q are a sample mean of 0.96 and a sample standard deviation 
©£ 0.21. The maximum likelihood method produces estimates of 
k with a sample mean of 3.93 and a sample standard deviation 
of 0.27. The values for q estimates are a sample mean of 1.00 
eeaea Sample standard deviation of 0.11. Although the method 
of moments and maximum likelihood method do not produce signi- 
ficantly different values for the mean of the estimates of the 
parameters, the lower estimated standard deviation for the 
likelihood method makes this technique more desirable from 
the standpoint of precision. No bias was evident in either 


estimation technique with the precision attained in the 


Simulations. 
The second case was the low correlation case with k = 4.0 
and q = 3.0. Here the distinction between the two estimation 


procedures is not as clear (Table II.E.3 and Figure II.E.2). 

The method of moments produced estimates for k with a sample 
mean Of 3.99 and sample standard deviation of 0.17. The esti- 
mates for q had a sample mean of 2.98 and sample standard 
deviation of 0.18. The maximum likelihood method produced 
estimates of k which had sample mean 4.00 and sample standard 
deviation 0.16. The q estimates had a sample mean of 2.99 

and sample standard deviation of 0.17. It is clear that neither 


method of parameter estimation enjoys a distinct advantage 


Oi 
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Results of Search for Maximum Likelihood Estimates 
mapa GLAR(]) Model> k = 4.0, ¢q = 3.0 


Value of 
Data Starting ee Run Time Number of Ending Value Maximm Difference 
Set k (in Min.) Iterations k gq Likelihood (1073) 
4.000 3.000 38 6 4.078 3.187 -620.637 : 
4.019 3.271 25 5 4.085 3.192 -620.637 
4.000 3.000 23 4 3.780 2.934 -655.488 : 
Pyvigy) 20543 26 4 30775 293d 6550459 
2 4.000 3.000 40 3 3.789 2.765 -663.886 : 
2 3.912 2.862 27 5 3.789 2.762 -663.886 
4.000 3.000 44 7 Woon aio o sesh 76 : 
3 4.000 3.116 50 4 Reals sill) | O55. 77 
4.000 3.000 44 4 ANG 0 eee 50" 78 
: 
3.958 3.174 68 6 4.152 3.385 -600.715 ae 
4.000 3.000 27 6 3.890 2.904 -638.667 
4.001 3.048 13 4 Bec snee oe Ge cnGcG ? 
4.000 3.000 43 3 4.051 2.873 -599.581 
4.221 2.902 49 E 4.049 2.878 -599.581 > 
4.000 3.000 23 2 3.954 2.934 627.062 
4.079 3.074 16 7 3.963 2.941 -627.058 Ee 
4.000 3.000 21 4 3.917 2.898 -637.234 
peraer 2.633 «32 6 3.908 2.891 -637.233 il 
9 4.000 3.000 15 2 4.111 2.924 -590.563 
9 4.093 2.906 43 2 4.121 2.934 -590.562 14 
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over the other with respect to precision or accuracy. How- 
ever, the method of moments is considerably cheaper with 
meeeect tO Computation time required. This is consistent 
femecmowh results that for i.i.d. Gamma data (k = q), the 
moment estimate for k is quite efficient when compared to 
the maximum likelihood estimator. 

paereeecase (kK = 0.75, q = 0.1825). The third type of data 
was a high correlation case with a low k value. Specifically 
me 0.7/5, gq = 0.1825, and the correlation was 0.75. As can 
mewseen in Figure II.E.3 and Table II.E.4, both the methods of 
parameter estimation considerably overestimated the parameter 
values, indicating considerable bias in the procedures. The 
method of moments produced estimates for k with sample mean 
of 0.8061 and sample standard deviation of 0.067. Those for 
q had a sample mean of 0.232 and sample standard deviation of 
0.026. The likelihood method produced estimates for k with a 
Sample mean of 0.852 and sample standard deviation 0.039. The 
corresponding statistics for q estimates are 0.224 and 0.004. 


As was true in the other high correlation case, the standard 


. Geviations of the maximum likelihood estimates are consider- 


ably smaller than those of the moment estimates. However, 
Since the evidence is that the estimates are highly biased, 
the advantage of this smaller standard deviation is not clear 
unless additional data would serve to reduce the apparent 
bias. The detailed results for this data type are presented 


meraote 11.6.4 and Figure II.E.3. It would be of interest 


94 





TABEE I.E. 4 


Results of Search for Maximum Likelihood Estimates 
in a GLAR(1) Model; k = 0.75, gq = 0.1875 


Value of 


Data Starting Value Run Time Number of Ending Value Maximm Difference 


Set k q (in Min.) Iterations k q  Likelihcod (1073) 
0.7500 Oe25 85 HOOl. = 02224. S25q005 
029795 0.2489 136 so65 1052229 3252655 
0.7500 Orle25 ag 652) “02235 5 cooeziG 
0.8026 Oa sisgh 19 Ose. §0-coo ss ooo.c 4 
0.7500 QeS25 12 . 806 .224 439.198 
Geez] 0.2147 66 Ate F/O PAR Ie ish e 
0.7500 0.1825 141 .866 0.222 1131.684 
0.7514 Ge22i3 5 2121 pOOO wm Oe222 LiSl. or 

. 7500 O25 Ss noOve BOe22o0cosoe Lee 
0.7749 0.2420 54 POO 12 ae On222 525 sean 
ee 00 Ope 2> Bo P90 0.217 eSaacou 

16256 Oeze20. 1210 POLOm eOe27 “Ces ecu 
0.7500 Oele255) 129 Seles) (OCP PR BUEN 
0.7308 O27 90 88 .884 UDO PASS SES, 
0.2500 O21.825 ya e795 me Os2235- “9772265 
OR 7759 Oez241, 117 195) 02223 977.265 
0.7500 eo25 66 .906 neo? 10182270 
0.8092 265 30 .906 P2227 10 eran 2 

. 7500 021325 ol 0.818 .212 410.842 
2535 0.2347 85 Qees2. 02225" 411 5703 
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to see if a technique such as (2-fold) jacknifing, such as 
that applied by Quenouille [Ref. 17] to correlation estimates, 
would help here. 

A much larger study would be needed to come to definite 
conclusions about the efficiency of maximum likelihood esti- 
mation in this model. However, as in the i.i.d. case for 
Gamma variates, m.l.e.'s are likely to be considerably better 
than moment estimations for values of k less than one. 

Note too that the maximum likelihood estimation did not 
include the mean value parameter. This could be done or the 
mean could be estimated from the sample mean X. The inflation 


of variation of X due to the correlation is known to be 





(asymptotically) 

a z 
Thus for pop = 0.75, the effective sample size for estimating 
u ina sample size of size n is n(l-o)/(lto). For 9 = 0.75 


Menis is n/7. 

The normality of the estimates was investigated with normal 
plots and Shapiro-Wilk tests for normality. A summary of 
results if given in Table II.E.5. The normality hypothesis 
is accepted at a 0.95 level if the Shapiro-Wilk statistic W 
is higher than 0.842, at a 0.99 level level if it is higher 
than 0.781. No strong indication of non-normality is indicated 


in any case. 
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III. MOVING AVERAGE MODELS 


fee LNT RODUCTION 

Although several researchers have proposed models for 
correlated, marginally exponential random sequences [Refs. 18, 
19,20|, Gaver and Lewis [Ref. 2] produced the first analytically 
and computationally tractable model for the generation of 
correlated, marginally Exponential random sequences. They 
showed that in the usual linear, additive, first-order, auto- 


regressive equation 
x = 8X + E (met Be) 


where {X, Melee 1S a Secona-order stationary, 
marginally Exponentially distributed sequence of random varia- 
bles; {E., n= 1,2,...} is an independent, identically dis- 
tributed (iid) sequence of innovative random variables; 

mee < 1, the distribution of the {E 3 which produces the 


desired marginal distribution for ice) is 


ie WlthepBobability 6, 

= = (Ci ie Aeee) 
- LE. iicialo) jeneeloyedesal jl abies, Vl jer 

where {EL, n= 1,2,...} is an iid sequence of Exponentially 


distributed random variables with the same parameter, i, as 


the {x} sequence. Equation III.A.1 can now be written as 
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(PXn-a with probability 6, 


an = GFriT.A.3) 
BX 1 + En with probability 1-8. 
haa Loa n= 1,2,...} is an iid binary sequence independent of 
a and {E,} such that ct 1) = Jos Tay(ih = a, eo firen 


equation III.A.3 can be more succinctly written as 


>, = BX + TIE. (iene A 4) 


x 1s a random linear combination of identically distributed 


random variables in the sense that the variable En actually 
enters into the sum only when the random variable I. has value 
one. Since for a given {I} sequence, the distribution of xX 


depends only on the distribution of x and En! ae will be 


ne 


Exponentially distributed whenever both x and En are inde- 


=k 
pendent and have the same Exponential distribution. This 
understanding allows the autoregressive relation in III.A.4 


to be transformed into a moving average by substituting another 


innovative random variable for the Xl to produce 


x =) pee + LA CE ilg 2 ) 


This model was designated the EMA(1) for Exponential moving 
average of order one. This EMA(1) model is one dependent in 


that xo and x are independent for j # tl. Consequently, 


+3 


only the lag one correlation, oes ie CORR(X_ ,% ), Or more 


n+l 


completely only the joint distribution of Xo and re need be 


studied. 
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Lawrance and Lewis [Ref. 5] give a fairly complete des- 
cription of this EMA(1) model. Of particular note is the 
Metative tractability of this’ model which enabled the authors 
to derive correlations, distributions of sums of X's, 
m@eensity function, spectrum of counts, joint density of 


x, and X conditional expectations, and other properties. 


@rodk’ 
The existence of these characteristics is beneficial in data 
analysis and is a primary advantage of the EMA(1) over pre- 
viously suggested models. However, the EMA(1) model does 
possess a limitation. The range of possible positive correla- 
=Lons , Py: 1s restricted to the interval from zero to one- 
quarter. Thus for a given correlation between zero and 
one-quarter, the structure of x and Xn+l and the sample path 
behavior of the sequence are determined. 

The structure of III.A.1 is that of a special random linear 
combination of Exponential random variables to given an Exponen- 
tial random variable. Other such random linear combinations 
are now known and for the first-order autoregressive case 
produce dramatic differences in sample path behavior of the 
- sequence {Xt. In this section of the thesis we investigate 
these random linear combinations in the context of the first- 
order moving average structure. 

In this Chapter we examine extensions of the model in four 
ways: 

i Negative Correlation 
McKenzie [Ref. 21] has suggested a scheme for inducing 


negative correlation in the EMA(1) process by correlating the 


16) ds 





sequence {I t. A better scheme, in that it produces a larger 
range of negative correlations, was introduced by Lawrance 
and Lewis [Ref. 8] for the extended first-order autoregressive 
model NEAR(1). This scheme involves a bivariate error se- 
quence {EES and its use is investigated in this thesis for 


the moving average process. 


2. New Exponential Moving Average Model of Order One, 
NEMA (1) 


It will be shown that no first-order moving average 
process which is a random linear combination of Exponential 
random variables can have correlation greater than one-quarter. 
Thus the differences in the processes for given correlation is 
investigated in terms of the joint structure of xo and X41: 
The NEMA(1) process obtained by using the NEAR(1) structure 
[Ref. 8] in a moving average context is analytically tractable 
and quantities such as the joint Laplace-Stiltjes transform 


of x and X Eae Spectrum Of Counts, P(X 7 4, and condi = 


gle aa aed 
tional expectations are obtained. It also combines the forward 
and backward EMA(1) models as extreme cases and is thus a 
natural extension of the EMA(1) model. 
3. The Moving Minimum Model 

A non-linear combination of Exponentials involving 
minima has been applied by Tavares [Ref. 22] to obtain a first- 
order autoregressive process which is intimately connected 
(Ref. 23] with the EAR(1) process of Gaver and Lewis [Ref. 24]. 


This structure is applied to the moving average process. It 


is shown that this process extends the range of attainable 


ILC 





correlations in first-order "moving average" process beyond 
that obtainable by random linear combinations of Exponentials. 
However, the process is slightly degenerate in that there is 
a positive probability that successive points will lie on the 


line Xe = bx, for positive correlations and on the curve 


-X =X 
+ , 
— +e — 1 for negative correlations. The price paid 


for the extended range of correlations is a limited analytical 
tractability as compared to the EMA(1) or NEMA(1) processes. 

The moving minimum process is investigated in terms of the 

fence Structure of XO and xX et peenough the 7OInNe aistribution 
Gan be derived, the functions are difficult to examine in 
detail. Therefore, simple characterizations of the joint 
Structure, in addition to the correlation, are examined. These 


include the P(X 22S I, a crude measure of the tendency of 


iho 
the sequence to exhibit runs up and down, and conditional 
expectations, E(X_|X,_,=x) and E(X__,|X,=y). 
4. The Beta-Exponential Model 

Finally, another random linear combination of Exponen- 
tials to produce correlated Exponentials is examined. Unlike 
» the previous models, the coefficients of the Exponential random 
variables are themselves continuous random variables. This in- 
creases the complexity and reduces the analytical tractability 


of this model. Simple sample path characteristics are derived 


or Simulated. These are special cases of the GIMA(1) from Chapter II. 


B. NEGATIVE CORRELATION IN MOVING AVERAGE MODELS 
The problem of non-negative correlations was addressed by 
meninzie [|Reft. 21} who modified the form of the EMA(1) model to be: 


OS 





Cli lB) 


where {E.. n= 0,1,2,...} is an iid sequence of Exponential 
fandom variables, {Zo, tow 15 a sequence Of binary 


fmeamcaom Vartltables with P(Z =1) = 1- P(zZ 


- aa) = 8, and Zo 


is independent of all En and all past Xo: By imposing an 
MA(1) correlation on the sequence {Zt, McKenzie was able to 
produce a negative correlation for the {xX}. However, this 
negative correlation is achieved at the cost of reducing the 
possible range of positive correlations for the {xX}. Using 


McKenzie's formulation, the range of correlations obtainable 


a as 
64'16°° 


An alternative procedure for producing negative correla- 


with the EMA(1) model is (- 


tions was introduced by Lawrance and Lewis [Ref. 3]. Their 
procedure requires two series of innovative factors {EL 
MeO, l....} and {El, n= 0,l1,...} which are correlated and 
may be antithetic. 

Antithetic variables are generated by using the fact that 


the variables E. and E defined as: 


(lela. Bile) 
where {U., i=1,2,...} is a sequence of uniform (0,1) random 


variables, are both marginally Exponentially distributed and 


correlated. P.A.P. Moran [Bef. 25] has determined that the 
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correlation between E. and Ey is approximately -0.6449 

ame this is the lowest correlation possible between 

Exponential random variables. E. and E! have a deterministic 
a 


relationship since 


(aie Bias) 


Using the Lawrance and Lewis extension of the EMA(1) 


process, the model becomes 


= ' 

x BE, + IEl a: (IIL.B.4) 
where {Ene n= 1,2,...} is an iid sequence of Exponential 
random variables, {Er n= 0,l,...} is a sequence of Exponen- 


tial random variables which are correlated with the respective 
variables in the Pee sequence, {I n= 1,2,...} is a sequence 

of independent binary random variables with P(I_= 0) Si = P(I_=1) 
oe or <p < 1, and Pease, {EB} are independent of each other 


and all previous x values. 


The correlation of the X's can then be computed as follows. 


Let E(X) = yp and recall that since Sa and {EL} have the 
Same marginal exponential distributions, VAR(X_) = VAR(E.). 
= t ' 
en (BE ait I yepP yp! (SE, t TLEL-1! 
Zz t QR t I q 
8 Be nh ene enn tn nel n=2 Snel nen n-1’ 
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Thus, using the independence of {EL}, cle and the iid nature 


of LPL and ae 


) = g“1.°+8 (1-8) [COV(E, ,B4) tu] +8 (1-8) w+ (1-8) “U7 


Zz | 
Ribs tock 8 (1-8) COV(E, ,E.) 


Therefore, 


COV (X14 /%,) = 8 (1-8) COV(E_,E) 
and using VAR(X) = VAR(E) 
CORR(X_ 4 /X,) = 8 (1-8) CORR(E_ ,E.) Chis B 5) 


Using Moran's result [Ref. 25] the correlation of antithetic 
Exponentials 1s known to be approximately -0.6449. Therefore, 
the greatest negative value that can be achieved for CORR(X_ 4 7%,) 
is approximately -0.1612 when 8 = 0.5. Since no restriction 

. was placed on CORR(E.,E!), the sequences {EF} and hee need 
not be antithetic, but can have any correlation that is possi- 
ble for two Exponential sequences with the same marginal dis- 
Semoucion. By specifying that EA = EL, the original EMA(1) 

is achieved and the correlation for the X's is §(1-8) as in 
Lawrance and Lewis [Ref. 5]. By allowing the correlation be- 
tween the {E_} and {E'} sequences to vary from -0.6449 to 


1.0, the correlation of the X's will vary from a minimum of 
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-0.1612 to a maximum of 0.25 (also depending on the value of 
8) as can be seen from III.B.5. The Lawrance and Lewis exten- 
sion of the EMA(1) model gives greater possible variation in 
the correlation than that of McKenzie, but requires two sSe- 
quences of Exponential random variables to achieve this range. 
Although it is clear that with EY = En a CORR(E_,E)) =] 
and when Ua) and ee are antithetic CORR(E_,E.) = -0.6449, 
it may not be obvious how to generate {Et} sequences with 
correlations between these two extreme values. A simple bi- 
variate exponential sequence with any desired correlation in 


the permissible range can be generated using the relationship 


with probability p, 
E Wievepiogadp lity l=p, 


where E; is the antithetic of E,. Then the extended EMA(1) 
model has two parameters, 8 and p, and the range of correla- 
tions for the X's is -0.1612 to 0.25, as above. The bivariate 
density for the pair {E;,E,} is not smooth. Other bivariate 
- densities such as thoSe in Gaver [Ref. 26] and Lawrance and 
Lewis [Ref. 27] can also be used. 

The above ideas on extending the correlation structure of 
the moving average models to encompass negative correlation 
can be applied to all of the new models given below. Details 


are not given. 


Oe 





fe Lee EXTENDED EMA(1) MODEL, NEMA(1) 
Po Lneroaucelon 
The original Exponential moving average process is 
discussed by Lawrance and Lewis [Ref. 5]. This paper 


considers the first order process defined by: 


| BE. Whehmomeobablility 8, 
SEL + Ent] Weweomoocdo li ity (l=3)., 
where {E_, n = 0,+1,+2,...} is an iid Exponential sequence and 


n 

eee. slis rendom linear Combination of Exponential 

variates is called the EMA(1) model for Exponential moving 
average of order one. Since x tsa LuneelLon. of En and Enel! 
this version is called the forward EMA(1). The backward 

version of EMA(1) is obtained when E+l is replaced by En-l 
Port .c.1. 

The fact that EMA(1l) is a single parameter model sug- 
gests that this model may not be sufficiently flexible to 
address all processes. Investigation of an alternate, two 
. Parameter model may indicate that a two parameter model is 
sufficiently more flexible to justify its increased complexity. 

The extended, two parameter model is based on the new 
Exponential autoregressive process of order one (NEAR(1)). 


The NEAR(1) model propounded by Lawrance and Lewis [Ref. 8] 


is defined as 
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| aa + BX ol Witimoneucdot lity & 
x = CEE Gari) 
c with probability (l-a) 
where {X, n= 1,2,...} is a second-order stationary sequence 


of Exponential random variables with parameter i, {E 3 is an 
imi quence OL Innovative Factors, 0 < 8 < 1, 0 <a< 1, and 
feo <« 2. By letting dy (Ss) and > _(s) denote the Laplace- 


Stieltjis transform of X and E respectively, Lawrance and Lewis 


eos r : : 
Secermaned that > (Ss) Be er 5S Closes iS arcane ater 1 cee 


fraction solution technique to invert this transform produced 


, oe oe 
| EL with probability TS (=a) 6 
7 (iE 23ers) 
f | Gio 8E When Deobabilit a RUC ea 
n Pp ye oes 
where VS n= 0,1,2,...} is an iid sequence of Exponentially 


distributed random variables which has the same parameter as 
the es sequence. 
By noticing that the autoregressive model given in 
Pewee ceustmig E1T-C.1.3 is a random sum of two iid Exponen- 
| tially distributed random variables, the NEMA(1) model is 
in the NEAR(1) model. 


meoducec by substituting E Oe aX 


n-1 
This procedure is identical to that used to produce the EMA(1) 


n-1 
model from the EAR(1) model and yields 


with probability ©, 
n ct (Cle bey. G)) 
E wieneonopability l=-a. 
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This model can be written in a more compact form as 
X = KE ee leis CRivG cc) 


where ae n= 1,2,...} is a second order stationary sequence 
of marginally Exponentially distributed random variables; 


{E n= 0,1,...} is an iid sequence of Exponential random 


13 ad 
variables with the same parameter as the {X_} sequence; 


(im n= 1,2,...} is an iid sequence of random variables with 


oud 
eee Poo — Oo) — a; {K, Hele) is anoiid se— 
quence of random variables with P(K, = 1) = 1-P(K, = (1-a)8) = 
oe) {I}, {K_} {E} lly j dent; 
ta) so? I, P KO , and En are mutually independent; 
fees lO < 6 < 1; and até < 2. 

The NEMA(1) model contains both the forward and back- 
ward versions of EMA(1) as special cases. When a = 1; 


P(I_=8) = 1, (l-a)8 = 0, and Engi nC) = 8. Hence, the 


NEMA(1) model becomes 


with probability 8, 


n-1l 
es | (lel Te Cr ey) 
SE W-1 + E, Wee WO TEoeaol lity )(l—6) < 
This 1s a form of the forward EMA(1). 
When 8 = 1; P(r =2) =a, (l-a)8 = (l-a), and 
= - = Beet Bee. = 3 
P(K, = (1 a) ) Toene 1. In this case, the NEMA(1) model 


becomes 


aehG) 





oo 


(l-a) EL Ween robabi lity (l=q) , 


X = (Gi ea a gs oon ee BY 


(l-a) EL + EL) Viemmomoanani Li ty oO % 


This is a form of the backward EMA(1) with 8 = l-a. There- 
fore, the NEMA(1) contains the special cases of EMA(1) when 
a and 8 assume specific values. 
One should also note in passing that the {xX} sequence 
becomes an lid sequence if a = 0 or 8 = Q. 
2. Correlation Structure 
The following relationships will prove of value in 


succeeding calculations 


P(I_ =8) =il- P(I. = 9) = a CEs Gee.) 
E(I)) =e (l—@je0 = ab. (ema G. 2. 2) 
P(K = 1) = 1-P(K, = (1-0) 8) Cie Gree) 
a ie 
> &(il-c)e- 
= ao8 Sy ee 
E(K.) = sel aiGe el a) 8 (3-7 y=ayR? 


1 Ba aC eee 
1-8+a8 


eee ae ac 6 - 
1-8+a8 


PaO eo 8 +cR ) 
I=6 a6 1-8+a8 


ae esl 





E(K,) = =o Cri i. C2 4) 


xX = K E +tH2E CT ies Gi) 
n neon non-l 
E(X) = E(K E+ IE 4) Saas 25) 
= E(K.E.) + E(I,E_,) 
= E(K_)E(E.) + E(I_)E(E_,) 
= (1-aB) E(E.) + aBE(E _4) 
E (X) = E(E) 
Since (os I and SEs are both Exponential, E(X) = E(E) implies 
VAR(X) = VAR(E). Since both {EL} and {x} have the same 


Exponential parameter, without loss of generality this param- 
eter will be considered to be l. This, of course, requires 
E(X) = 1 and VAR(X) = 1. 

The possible range of correlations for the NEMA(1) 


model can be determined by a simple calculation. We have 


Seal Antinne? * In4i®n>- 


Thus 


eZ 





wae 
~ 
| 


(KE ste liegt ) (K 


n n+l n n n-l etre n+l on: 


K + K mao log Ae E 


eee on eee ten 1 


ne EE 


eon nn-l 


Using the independence of {kK}, nee) and {E} ancmene 91d 


nature of these three sequences, we have 


E(X.X,.,) = (1-a8)*(E(E,)1* + (1-a8)a8 (E(B) 1° 
+ (1-a8)a8(2{E(E_)}7] + (a8)? (B(E,) 1° 
=e (=a) as 
Therefore, 
COV(X_,X_ 4) = (1l-a8) a8 
and 
CORR(X+X_,,) = (1-a8)a8. (ea 


Therefore, the original NEMA(1) model has the same range 
of possible correlations as the EMA(1), namely that the corre- 
ae One should note that it 


is not possible to distinguish the parameters %@ and & from the 


Macions must lie in the interval [0 


correlation, or even whether the product, a8, has a given value 


Or one minus that value. This iS similar to the normal tmnoving 


3 





average model of order one where the cases 9 and 1/6 are 
indistinguishable. In the normal case, the range of 9 is 
limited to the interval [0,1] on the basis of invertibility 
[Ref. 28]. It would seem simple and convenient here to limit 
a8 to the interval (0,51. However, non-normal processes are 
not completely determined by their correlation structure. In 
fact, Jacobs and Lewis [Ref. 6] showed that in the EMA(1) 
case, the values 8 and (1-8) can be distinguished using direc- 
tional moments, EMXeacte and BQO ae Hence, such a re- 
Striction on the value of a8 is inappropriate. 

One should also note that, since the correlation 
between x and X 4K is zero for all K with absolute value 
greater than one, the first-order correlation completely 
determines the correlation structure of the model. 

The restriction on the range of attainable correlation 
is disappointing but not surprising since it can be proven that 
any Exponential moving average process of order one generated 
as a linear combination of independent Exponentials must have 


a correlation that lies in the range 1 pedie The proof of 


-this contention follows the previous calculation closely. 


fo SOREM: 

Assume oy We or eres a Sequence of iid Exponen- 
meal VYariaoles with unit mean, {A,, = lees. oy aril {BL 
n= 1,2,...} are sequences of iid random variables, and {AS, 


rt, {E 3 are all mutually independent. Moreover, assume the 


sequence eae n= 1,2,...} defined by 
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Xx = AE oe Jo) 6 LT £2C 216) 


is a unit mean Fxponential sequence. Now 


E(x) E(ALEA + BE 


n oe? 


E(A,)E(E,) + E(B,)E(E,_,) 


- 
il 


E(A,) + E(B.) Ci 2 a) 


imeaddition, since Xn PO E@useal lei. both oon and BL, must be 


non-negative for all n. Hence E(A_) > 0 and E(B.) 2 ae 
also follows that l > E(A) and 1 > E(B). Now 
eel ~ (AE t Bene nar) (Aner Png + Bnd En! 


Therefore 

_ 2 
ee n+l! 7 nent lon areienendian=1° “none no 

% Bs een n- 1) 
2 2 

= [E(A)]~ + E(A)E(B) + 2E(A) E(B) + [E(B) ] 

= [E(A) +E(B)]* + E(A) E(B) 
Since E(A) + E(B) = 1 from III.C.2.9, E(A) + E(B) = l. 


eee CA) E(B) 


aN Es. 





Therefore, since E(x) is one, by assumption, 


COV (XX, ) = E(A)E(B) 


+] 


E(A) [1 - E(A) ] 
and 


CORR(X_,X, E(A) {1 - E(A) ] (ian G2 eee) 


+1) 
BemeeeO = 4A) < 1, the correlation must lie in the interval 
[0,5]. Q.E.D. 

The possible range of correlations can be extended by 
reformating the model. We choose to do this first by using 
the method devised by Lawrance and Lewis [Ref. 8] and given 
in equation III.B.4. With variables and sequences defined as 


mamecwation III.C.1.5 let {E} and {Et} be correlated sequences 


and redefine the NEMA(1) as 

me el ee ee RR (Erie C222 ae) 
non n 

Then it follows that 


= ’ 
ax Cree L E ) (K 


no n+l pee ne? © on+1 ow 


= &«K K = raat = ee Pee 


q 
en ers ne lon ntl n-1 


{ | 
as eon a ne 


IgE 





Thus 


E(X X41) = (1-a8)? + (1-a8)a8 + (1-08) a8 (COVIE,,B']+1] + 078° 
= et (1-a8)a8COV(E_,E!) 
and 
COV(X. X14) = (1l-a8)a8COV(E. ,E? ) 
Fina ly 
CORR(X_, X14) ~ (l-a8)a8CORR(E_,E,). CEIMmG 52.12) 


As described above, Moran [Ref. 25] has determined that 
the range of possible correlations for two Exponentials is 
(-0.6449,1.0). Thus whenag = 0.5, the possible correlations 


Zeng x and xX Puede tie anterval (-0.1612,0.25). This 


aoe 
procedure extends the range of possible correlations at the 
cost of generating the additional 1BS J sequence. 

McKenzie [Ref. 21] has suggested that the range of 
correlations could be extended by requiring that the ce 
sequence be correlated. Using this scheme, he was able to 


show that the correlations for the oe sequence lies in the 


interval (- grag) Because of the requirement of the moving 
average process of order one to have zero correlation for lags 


greater than one, the correlation of the iri sequence also 


Sha ae 





Maemtemoe OF MA(I) type. It is this restriction that pro- 
duces such a narrow range of correlations. A logical and 
obvious extension to McKenzie's work is to reguire that both 
the Ke and ore sequences in the NEMA(1) model have a MA(1) 
correlation structure. This can be combined with the corre- 
lated {ELI, 1ED scheme of Lawrance and Lewis. If this com- 


bined case 1S carried out, the NEMA(1) model is as follows 


ns Gets AT te (it G~ 2eieot 
n non non-l 
where {KL, n= 1,2,...} is a sequence of random variables 
with an MA(1) correlation structure with Pk = ©) = 
» ie _ a 1-8 : a 
1-P(K, = (1-0)8) = y=7i=g7R! {I,, n= 1,2,...} is a sequence 


of random variables with an MA(1) correlation structure with 
mon. = 6) = Lae ot 0) = a; oe and {Et} are correlated se- 
quences with marginal Exponential distributions with unit 


means; and (KI, UE er and {EL} are mutually independent. 


Now 
= t ' 
Ponon+1 (Ket I yEn-1?) Ryer Ens * Ing1®y? 
= ' ' 
een Se ten nelene ll Sn ns] nn 
] t 
- en ae na] 
SO 


IES 








n= BECK )+(l-a8)a8+(l-oB)aSE(E E!)+E(I_ I.) 


aeoless 


n+l +1 
i — 2 = 
= COV(K 14 /K,)+(1 eo) +(1l=a68) ab 
; 2 
+(l-a8)ag Sey COV AL Ete) 
= = ! 
1+COV(K, 4, -K,) +CoV(T iT) +1 aB)aBCOV(E. ,E) 
Therefore 
COV(X, X14) = COV(K, 5 ,/K,) +COV(I. 4 I) +(1-08) aBCOV(E, , ES) 
and 
CORR(X_,X__4) = COV(K, 4, ,/K,)+Cov(I_.4 I.) (ae? 2a) 


+(1-a8)aB8COV(E_,E) 


Although this scheme obviously extends the range of possible 
correlations, it does so at the expense of considerable com- 
plexity. Considering the limited range of correlations 
possible by imposing a correlation on the ae and OSE 
sequences, the additional complexity may not be warranted. 

If in spite of the complexities involved, one decides to in- 
duce correlations in the coefficient sequences, the NEMA(1) 
because it has two such sequences will yield a slightly larger 


range than the EMA(1) model. 


JRE) 





One of the possible advantages of a two parameter 


model is the capacity for modifying P(X >X_) and, conse- 


leet 
quently, the sample path behavior of the process while main- 
taining a constant correlation. Since the correlation is a 
function of a8, one can vary the values of a and 8 while 
keeping the product constant. The P(X 44,7 X,) can be calcu- 
lated by addressing each of the sixteen possible combinations 
of K and I values for X41 and Xe computing the probability 
for each combination, and weighting the probability associated 
with a given combination by the probability that the given 
combination occurs. A sample calculation is provided and 


complete results presented in Table III.C.3.1. 


Example: We have 


<> AE us eaene Tt! 
eros ntl ntl ~ on+i”n’ 
ieee I=P(K = (1-c)8) = ++ 
n n | l-(l-a) 8’ 
P(I_= 8) = te a= OY ea! 
Let + = 5 = 6 and consider the case where 
eo yet = 8, K, = 21 





Since the ie and Glee sequences are both iid and 
independent of each other, the probability of this combination 
of parameter values is simply the product of the individual 


probabilities of occurrence. Hence the probability of 


@ecurrence is oe Then in this case 
Bs aa, ead 7. ree ye 1! 
= P(E 4, > (1-8) EL t+8E__y) Cl ae Oat ee 
Now E_,, 1S independent of Y = (1-8)E_+8E__,. Therefore, the 


calculation required by equation III.C.3.1 is straightforward 








once the density of Y is obtained. We have, with fe (x) 
ghar dl 
Bucep.d.f. Of FE, (i.e. e -) 
y/8 
P({1-8)E,+8E,,<y) = J EUIATE Ie SVE, “fp (x) ax 
y/8 
= f Gg Any sors B = Xx) (x) dx 
0 n-1 
y/8 
= Via x 2 
= - P(E LS T= [Fyn =) fg (x) dx 
y/8 So ue: 
Se eteete: 8 )e. “dx 
0 
Bx 
y/8 = y/8 ay [x-q=~] 
=f e dx - f ee Im8 "ax 
0 0 
ne Sais 
z Se S/S ese 
Sree i e ax 
0 





P([1-S]E +8E__j <y) 











/8 TF 2 ter 
s _ ‘a = 
l-e es i j” ob 
~¥ ;in28 
= eae 8° 1-8 
l-e said Oa 5] (1- = 
/8 1-8 “5 te le 
Se ee : = ae 
- Serge oe Sieg 
8 -y/8 1-8 ~T=8 
ie 4 < as Tragli-e I ae 1X (dle ) 


Therefore, the density function of Y, fy ly), is 


{a -y/8,,1-8,, 1  -y/(1-8) i 
(tq) (e)e + (Fooq) (PS = , for 8 # 5. Gaver and 

Lewis [Ref. 2] gave the necessary and sufficient conditions 

~A,% “A5x 
for a mixed exponential of the form T,A,¢ +7,r,€ ; 
T) ae T4 + T 5 = 1, and Ay < Ao to be a proper density 
Memctlon. The condition is that 

lel 
aT = (Ali = CLG eg rece 
i- Ay 


In this situation we address two separate cases. The first 


Is 
i=2s 





case 1s when 0 < 8 < i. In this case > 1 and the 


2 


requirement is 


ie oe _ 8 etl 
Toe 1 tt = Ca) rie) oa ee a 
ie 
= IE2E 


ae 








ene LEE.C.3.2 is satisfied. 


ime tenis case - — > 1 and the requirement is 
WAL Ah) ee 
-8 a 
lel: B J 
aye: 

=e ea 
fea TTtT.C.3.2 1s satisfied. If 6 = Sr the density of Y is a 
Gamma(2) density. Therefore, the p.d.f. of Y is a proper 
density. This result can be used to complete the calculation 
Gq P(X417 %,)- Recall that equation IITIC.3.1 stated 

coe a? 7 eae)” ee 1 
= P(E 417 ¥) 


my conditioning on 





Varcdmedeusing che p.d.f. 


The second case is when i — 6 <s. 


2 


for Y derived above 


this 1s 
P(X.1>%,) = J PEns1 > S797) a iehy 
- V 
= -y/8,, 1-8 1 1-8 
= J e- (7 ioe) =)e Gree eae © J dy 
B+1, 
-_—, 2 (Fe S 
= eeyey gens Sot 
oe) (228) 
ile 1-8 2-8 1-8 
+ (759) (Zag)! (7=R)e dy 


es 
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Table I1I.C.3.1 presents the results of the calculations for 
all of the sixteen combinations of parameter values for Xe] 
and x: 

When the P(X 44 7 X,) was calculated for various values 
Soa ana 8, Lt was found that the values for this probability 
varied from 0.44 to 0.54. Table III.C.3.2 contains the results 
of these calculations for four hundred forty-one combinations 
of parameter values. Although the variation in probability is 
not large, it does represent an increase over the forward 
EMA(1) model. In particular, the forward EMA(1) model can not 
produce a probability greater than 0.50. Consequently, the 
NEMA(1) model not only has a greater range of possible proba- 
bilities, but also can produce probabilities greater than 0.50. 
The implications of this greater range is that the NEMA(1) model 
can address data sample paths that have a slight tendency for 
 @lither runs of increasing or decreasing values, while the 
EMA(1) can only address sample paths that tend to produce runs 
of decreasing values. 

Examples of scatter plots and sample paths for three 
sets of parameter values and positive correlations are given 
in Figures III.C.3.1-III.Cc.3.6. Because of the relatively low 
correlations possible and because of the limited range of 


values for P(X 4,7 X,) differences among the figures are not 
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Sharply delineated. However, differences can be detected, 
Bazticularly in the sample paths. In Figure III.C.3.1 with 


an @ value of 0.95 and 8 value of 0.50 has a P(X a) = 


pho dl 
0.44, the lowest value for this probability. Since this pro- 
duces a slight tendency for runs of decreasing value, the 
number of extreme values (i.e. greater than 3.0) 1s two. 
Boeraaures IlIT.€.3.2 and IIIC.3.3 the NOSE ue Se is 0.50 and 
0.54, respectively, with a corresponding increase in 
the number of large values. This trend is more 
Micticule to detect in the corresponding scatter plots. 
Piaqures III.C.3./7-I1I1.C.3.12 provide sample paths and scatter 
plots for the same a and 8 values as previously displayed, 
but with antithetic innovative sequences (see III.B) and 
consequent negative correlations. Although the negative 
correlation is evident, trends in these figures are difficult 
to detect. The extremes of sample path variability produced 
by the NEAR(1) process [Ref. 8] are not reproducible with the 
NEMA(1) process. This may be attributable to the restricted 
range of possible correlations. 
Seba tace Transform of Sums 

One aspect of the EMA(1) model is its analytical 
tractability. This is evidenced by the ability to derive the 
Laplace transform of sums by a recursive relationship given 
by Lawrance and Lewis [Ref. 5]. This tractability carries 
over to the NEMA(1) process. The Laplace transform is useful 


in obtaining quantities which are of use in analyzing point 


processes, namely the intensity function and the spectrum of 
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counts. These quantities are derived for the NEMA(1) process 


in subsequent sections from the results obtained here. 


ieee ederinedvas Gi irr.c.?.5, let —=—-—. = § 
n - 1-(l-a)8 ist 
and (l-a)8 = y. Further, let } X, = T, and ¢, (s) = E(e “). 
i=l ie 
Then we have 
cE = Xy - Xo + ... + x oe ie ey Pe aS 
= KE) + T,E, > K5E. + Io, ae igo KUL + TE oy 
Sere te KE e+ le) = eee 
Then letting L. = I. TeV ele oro ana using Ehe 
J eral J 
mutual independence of the iid sequences eee DORE {EL} 
-sT, 
Om is) = E(e ) (iia Gesu) 
xr 
: ees 15 1 -+L,E,+I,E)] 
-SK E sL E -SL,E -SI,E 
aera 16 TB (e r-l'r 1) E(e Ii 1) Be IL oe 
-SK.E. —-SI.E. -SL.E. 
Now let Ys) = E(e J Je #7 (s) = E(e J J), #7, (Ss) = E(e J Jy, 
Then 
Pe aye (s)tre is) 1° Gti Oa 3 
TY K i > IE ° ° 


To evaluate these guantities note that 
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IL Witmemorobpability 6, 





K = 
y When OnObapalaty l=. 
Then 
Se 
¥y (s) = E(e ) 
=—SP,. SNE 
= §E(e J) + (1-65)E(e am 
¥(s) = 56, (s) a (1-5) o, (Ys), 
where $¢.(s) = — So 
E l+s 
= 6 (1-6) 
¥x tS) ae aaa 1l+ys 
Also 
( 8 With orobabllity a, 
I = 
0 Wen Probability l<a, 
so that 
-sI.E.) -sRE 
¥i(s) = Efe JJ = gE(e Be (ea) 
¥ — a —_ 
7s) ~ L+e8s + (1-a) 
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(a A ea Oa ey 


(Cl Da Ge 4ec6 


(ere cee) 





To evaluate ¥, (s) note that 


B+1 with probability ad, 


B+y with probability a(l-6), 


iy = 
Ei Wien probability (l=a)<d, 
Y with probability (l-a) (1-6). 
Therefore 
~sL,E.) 
¥7 0s) = E(e 
-s{8+lL]JE. -s{8+vy]JE. 
= adE(e J) +a (1-6) E(e J 
-SE. = Se 
mma ror te. -) + (l=a) (1-6) Ele J) 
is) = a6¢.({8+1]s) +a(1-5)o,([8+y]s) + Gio os) ert te.4 8) 


+ (1-a) (1-6) 6, (ys) . 


Meeemen tne results of I[I1I.C.4.3, I11.C.4.4, ITII.Cc.4.7, and 


mr. C.4.8 


dm (S) = [69,(s)+(1-6) 6, (ys) 1x lad, (8s) +(1-a) ] (a a 
xfado,([8+1)s)+a(1-d)o,([8ty]s)+(1-a) oo, (s) 


+ (1-a) (1-6) 6, (ys) 17+ 
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This extends the result 


to the NEMA(1) process. 


Bo. Eaerace Transform of the Distribution 


SzVFCounts 


eee ay 


in Lawrance and Lewis 


Pret. 5] 


The Laplace transform of the sum is useful in deriving 


the distribution of the synchronous counting process of the 


number of events that occur in (0,t] when the origin is estab- 


lished at the occurrence of an arbitrary event. 


of events in (0,t] 


the relationship 


ie 
where Ny 


of the first interevent times. 


P(N, =r) 


where P(t) Homeieomo Seri bution function of TL. 


bility generating function of we 


Thus 


t 
me 
¥¢(Z;t) =e (02 =) 
_ 
> a Eas — ©) 
alse 
= ate Ut Ge} F 41 (t) 1 
= r-l 
= 1+ (Z-1) ) 2 F(t). 
r=1 


The number 


is related to the distribution of a sum by 


(ai C2 oe) 


is the number of events in (0,t] and T is the sum 


The proba- 


can then be written as 





* * 
Let ¥ ¢ (278) be the Laplace transform of ¥e(Zit), and f(t) 


the Laplace transform of f(t), ac O.d.f. OF ee Then 





¥. (238) = = . &1=2) d oe (III.C.5.2) 


Using the Laplace transform of the sum from III.C.4.9 





(1-2) 5 gr-t 


l 
CO 
l 


* 
¥ (278) [9d (s)+(1-5) oO, (ys) ] 


r=l1 


x lag, (8S) +(1-a) 1x (06d, ([8+1]s) +a (1-8) >, ([B+Y 1s) 


+ (1-0) 66,(s)+(1-a) (1-6) o, (ys) 1777 


¥. (278) = =~ $422) 15 4.(s) + (1-8) op (ys) 1* lad, (85) (ance lene 
LMS) 
1-z[ado,([8+1l]sta(1-6)o,([8+y]s) 
+ (1-0) €4,(s)+(1l-a) (1-5) ¢, (ys) ] 
where o_(s) = lt 
E l+s 


ce me (t) fee mimEenisl ty funetion of the point proc- 
. ess, then me (t) , its Laplace transform, can be obtained by 
differentiating III.cC.5.3 with respect to z, evaluating the 
derivative at z = 1, and then differentiating with respect to 
Ss. These steps, when taken, produce a series of tedious 
calculations which produce no analytical insights. The result 


of these steps is 
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eerie ona (se453s=20b-—Sa8-4+0-8~-) s- 


MT Os 508 --30e> +o -B-+a-8-) s° 
een ea ee CS PF II.C.5-4) 
s+(1+48-3a8+a' 8° )s°+(38+58° -2a08-70a8 +3a°8°)s 


+ (28°+28°-308°-308-+a-B-+a-B-)s" 


This result can be verified in a number of ways. First, when 

a = 1, the process is the EMA(1) process and, hence, the formula 
must reduce to (4.2) given in Lawrance and Lewis [Ref. 5] with 

XK = 1. Second, when a = 0, the NEMA(1) process reduces to a 
Poisson process and the formula under this condition must 

reduce to the Laplace transform of the constant intensity 


function of a Poisson process with rate l, =. aa py ew teh 


8B = 0 the NEMA(1) process iS again a Poisson process. Finally, 
* 
using one of the Tauberian Theorems, lim m-(t) = lim sm, (s) = 1. 
t+o s>0 
We take these cases in turn. First, when a = l 


Tena = one 52 -—raeeisue-sa-ec)s- 


a - +(28°+28°-308--308-+a-B-+a°B-)s> 
me (S) a 


5+ (1448-308+0-°8~) s-+ (38+58 --208-708-+30°8~) s- 


PE ee, = sa- 308 Fa°k-40-B") s- 


reduces to 


1+ (1428) s+(8+8°)s* 


* 
m-(s) 
f Reivers) s-=(e+8-)s~ 
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* (1+8s) (1+[1+8]s) 
m.(s) ——— 
s([l+s] (1+ (6+87)s]) 

(1+8s) (1+[1+8]s) 
Eeclos (1 +s) ( 


? 


a 
a lewayy | 3 


which is the result of Lawrance and Lewis [Ref. 5] with A = l. 


In the second case with a = 0 


iPaieeaianetn (3456 -ha-+(26°+28-)s~ 
5+ (1+48) s~+(38+58°) s-+(28°+28-)s" 


* 
m- (Ss) 


i 
S 


? 


the Laplace transform of a Poisson process with rate of l. 


In the third case 8 = 0, so 





gaan the Laplace transform of a Poisson process with a rate 
oc 1. 


In the final case apply the Tauberian Theorem 


Pa pbeeeniayem(sess5e--20e-5ae-+0-6-)s- 


; * + (28°4+28°-3a8-+0-8~+0°B°)s> 
ee ie ee 
s>0Q 1+ (1+48-308+a 8 )s+(38+58 -2a8-7a8 +3a 8° )s 


+ (28~°+28°-3a8-+a-B-+taB-)s> 


HH 
= 


as required. 
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6. The Spectrum Of Counts 


For the statistical analysis of series of events the 
most useful quantity associated with a process is the 
(Bartlett) spectrum of counts. The spectrum of counts, 
G,(w), 1s the Fourier transform of the covariance density of 
N-(t). It 1s related to the Laplace transform of the inten- 
SEty function, me(s), by the relationship derived by Cox and 
Lewis [Ref. 29] 


_ x cane aoe 
g, (w) = ~(Ll+m,_fiw] +m,[-iw]). 


We now derive this for the NEMA(1) process using III.C.5.4. 


* 
In that expression for mo(s), let 


a, = 1 + 48-208, Chink ee 

a5 = 38 + 58° - 2a8- See Ame ae, Cie. 

a, = 28° ~ 28° = 308° - ee - 0°3* + ace ee ee ae 

by = J] + 48 - 308 + rae eumaren 

bo = 38 + 58° - 2a8 - Wee + Ree BS, (Che ee 

b. = 28° = 22° = ae - Bice + “eae + sees (Thies Gr 
Then 


a7 


ly) 


2 


4) 


5) 


6) 





A lta,sta,s“+a,s> 
nay a SR 


2D 
s+b,s +b.s +b.s 


Recall that A = 1 so 


i lta, (iw) +a, (iw) “+a, (iw) > 


Aire 


= (yw) = oor 7 
ui iwtb, (iw) “+b, (iw) *+b, (iw) * 


lta, (-iw) +a, (-iw) “+a, (-iw) ° 
-iw+b, (-iw) *+b, (iw) +b, (-iw) * 
[iw+b, (iw) “+b, (iw) ?+b, (iw) 7] 
: : 2 3 : 4 
x{-intb., (-iw) “+b. (-iw) ~+b. (-iw) *] 
a. a Ee ee Chil .Ci6na5) 


g,w) == 
“ i [iwtb, (iw) “+b, (dw) ?+b, (iw) 71 


« [-dw+tb, (-iw) *+b, (-iw) ?+b, (iw) 7] 


[lta, (iw) +a, (iw) “ta, (iw) 7] 
x [-iw+b, (-iw) “+b, (-iw) *+b, (-iw) 7] 


[iwtb, (iw) “+b, (iw) ?+b, (iw) “1 


« [-iwtb, (-iw) ?+b, (-iw) 24b, (iw) “] 


[Lta, (iw) ta, (-iw) “ta, (-iw) 7] 


x [iwtb, (iw) “+b, (iw) +b, (iw) 7] 


[iwtb, (iw) *+b, (iw) “+b, (iw) 7] 


«{-dwtb, (-iw) “+b, (-iw) *+b, (iw) 7] 
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Consider the first term with numerator and denominator the 


same. Let the denominator = D. 


D = [iw+tb, (iw) “+b, (iw) +b, (iw) *] [-iw+b, (-iw) ?+b, (-iw) ?+b, (- iw)“ 
= w°-ib w>-bow*+ibw°+ibw>+b{w ib )b,w°-b bw bow" 
+ib bow? +b5w°-ibsb4w/-ibjw°-bjb4w°+ib,b,w +b4w° 
= i (-b,w> +b jw? +b, w>-b, bw +b, b,w?-byb4w/+b5b,w /-b aw”) 
+ (Ww “bow +b Sw -b) bw -bow +b4w°-b, baw +bsw®) 
= w*[1+[b%-2b, Jw°+[b3-2b,b,]w +tbSw°) (III.C.6.8) 
Let 
X = 14(b4-2b.)w*+(b$-2b)b,)w tb6u, 
where bo, bo, and b. Paeemderineduay  011.6.6,4, ILl.€.6.5, 
rir .c.6.6, respectively. Then 
DS Meer e668) 


Consider the numerator of the second term in III.C.6.7 and 


@all it N2. Then 
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a 2 eee ae 2 - 3 ue 4 
N2 = [l+a, (iw) +a, (iw) +a, (1) ll iwtb, ( iw) +b. ( ao +b, ( Teh | 
ee Die 3 4 Ze Sis. Ab se 3 
= 1W bjw +1bow +bu ta,u 1a,b jw a bow tia ,b.w tla.w 
4. a 6_ de. 5 6. 7 
tab) 1lab.ow abu au tla3zbjw ta .bow 1a,b.W 
N2 = i(-wt+la.-a.b,+b Neha b.-a.b.~ta.~b ene b ts (lebiee. 6 2.10) 
eee = 2 eer pee? os). B73 Et 


- 2 2 _ 4 _ 6 
+(ay bj )w +(a5b, a,b, a,tb,)w + (a,b, aob,)u ; 


where Ayr Apr az, b b anGeoovarereertined 1n LTIt;C.6.1 


ged 3 


through III.C.6.6 respectively. Consider the numerator of 


mie wcoLra tcemm in Iff.C.6.7 and call it N3. Then 


i — s = e 2 — * 3 LJ e 2 8 3 8 4 
o6|CU= [l+a, ( iw) +a, ( avs) +a, ( iw ) ] [iw+b, (iw) tb. (iw) +b, (iw) ] 
a. 20m ae S 4 eee 3 4. oa 3 
= 1) bjw 1b6w +b3u ta,w t1a,D)u a bow 1a,5D.0 1a5u 
+a.b ears b oe b es res b Mees b home b i 
2° i De Dy 3s S Bicol: See Sages 
ee ee ea pen (ae =a.b.+a.b,lw°+ta-b.w’)  (III.C.6.11) 
7 ee as lag? ee 22 3" At a3 oo 
+(a,-b yw+(a baeaebe- ast oa Deca ae 
gs Fe Zoli Suns: Ca “263 : 


where Ayr Agr 43, Die Do, and b. are defined in III.C.6.1 
Ehprough T11.C.6.6, respectively. Note from III.C.6.7 that 
all terms in the sum have the same denominator. Use III.C.6.10 


and III.C.6.11 to determine the numerator of the second and 
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Berd terms and call it N. Then 


7 a ” os =o | vig ~ 
m = if w+ (as a,b ,+b,)w +(a,b alia. +a b,)w a,b Ww +W (a, ajb,+b,)w 


~ (a, 53-a5b, 


+(a,b5-a5b, 


Se b,)w +(ajb,-a b a,tb,)w + (a,b, ajb,)w i” 


i 


het 


Z 
y = (a,-b )+(a5b -a,b -a,tb,)w +(a,b,-a,b,)u iP. 


saa 


Then 


Meng LlL.C.6.7, 


g, (w) 


where x and y are 


respectively. 


-a.b.ta.b 


3 


Se a: 3S 


5 7 2 
32,4 tazbiw ]+(a,-b,)w +(ajby-a,b, 


6 I | 4 
Yw +(a,-b,)w +(a5b,-a,b a,7b,)W + (a,b,-a5b,)u 


2 


6 
sae 


6 
prig ae 2 


ihe. Ono ott. C.,6.13 





Ly xX a 2W) Y) 
{oe 2 

i) WwW xX 
1 (Xt2y) 
TT x ‘ 


Gominecuimme Eieec.o,. 9 and ITII.C.6.12, 


Beds 


4 


-a.t+tb.)w 


Bane: 


6 


(Pie. 6.2) 


pill a Bos ra ci ee) 


Crime 46 i414) 
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Figures III.C.6.1 through III.C.6.3 show the results 
Seeene Calculation of the (Bartlett) spectrum of counts. In 
presenting the results the constant = in III.C.6.14 was ig- 
nored. Figure III.C.6.1 shows the spectrum of counts for the 
Same a and 8 values that were used for the sample paths and 
emerer lors Of Figures T11.C.3.1 through IIIC.3.6. This 
figure also shows the variation in the spectrum of counts as 


the P(X, > x) varies from its lowest to highest values. 


+] 


moire T11.C.3.2 holds the P(X, » Xx) constant and varies 


a 
meme correlation. Since the spectrum of counts for a Poisson 
process 1S a constant one when A equals 1 and the constant 

= is ignored, the correlation can be viewed as a measure of 
the process' departure from a Poisson process. This 
divergence as a function of the correlation shows clearly in 
this figure. Figure III.C.6.3 holds the correlation constant 
and varies the oS sesh Re The slight variation in the 
spectra shows that while the spectrum of counts does vary 


with the P(X, > Xi) the correlation plays a more dominant 


+] 
BOLE. 

The analysis from the Laplace transform of sums in 
III.C.4, through the Laplace transform of the intensity func- 
Meme nett? .C.5, =o the spectrum of counts in this section 
can be performed using the correlated {Et {z,} sequences 


of III.C.2 and thus for negative correlations. Details are 


not given. 
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Peeve ee eaplace-Stieltijes Transtorm of x, _and X oe 


Because the NEMA(1) process is, by construction, 
only one-dependent, all of the second-order properties of 


{x} are contained in adjacent pairs (OS pee a ee In previous 


sections quantifiers of the distribution of {X_ 1X Pesien 


n+l 


as oF and P(X > xX) have been derived. Here we give the 


‘ole ald 
Laplace-Stieltjes transform of the joint distribution. One 
could, for example, study the effect of the two parameters 
from this result by deriving directional moments. 

The joint Laplace-Stieltjes transform of Xn and Xt] 


can be calculated by considering each of the sixteen possible 


combinations of parameter values for xX and Xe]! as was done 


a 1-8 xs aaa — a — 

mr bit .C.3. Let Lane = On, (l-a)8 = { Py es (Sj So) 
“S41 X,7S80%n41 l Ween ek 

E(e A ebars! op (s) = Then 


> 6 





aoe nee, |-SalEe {toe ) 
6 (eee. ) acne pr n-l 2° ntl n 
oe, X ee 2 
mw n+l 
Soo eo Beal Se 
+ (l-a)assE(e + 2 B-l 2 ntl) 


-s,E_-s.[E +BE_ | 
Meiewisce(e o 2 4 ntl 'n 


~S,E -S.wAE 
Gemmdeaicete 6 2 * ntl, 


+ 


-s,{[yE_+8E ]-s[E +BE_ ] 
ase eS ume 1 n n-l 2 eri | n ; 


-s, [YE_+BE ]-S.E 
+ (l-a)adé (1-6) E(e 1 a n=t ye Dk 


-s,yYE -s.~[E +BE ] 
MGs (l=s)Ete + "™ 4 nti ny, 


“~S,YE_~-S.~E 
PG Orieostl=s\e(e © 2 * mt, 


-s,{E_+B8E ]-s,[vYE +BE_ ] 
+ aa(1l-é6)dE(e Ien n-1 2 avs . ) 


poeta to toe | say E 
Pee iiiesjisece © & n-+t 2 nti, 


Bere oa (Vy bakt3E. | 
+t a(l-a) (1-8)sE(e 22 2 nti nn, 


cas 


Si n 27" n+1, 


Pero mela e( l—o) 6B (e 


See ee ees Ly 2 +RBE ] 


Sou ( noobs . |=S¥E 
Pees i-sjE(eo + oD Nl 2 ntl) 


-S,YE_-s5[VE,_,+8E,] 


+ a(l-a) (1-6) (1-S)E(e ) 


—-S,VYE_“S.AYE 
+ (1-a) (l-a) (1-5) (1-8)E(e 2 7 2 MI, 
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(S/S) = aad6o,($,+8S,)€,(8s,)€,(s,) 
(1-a)ado0o,(s)) o,(8S,) >, (s,) 
a(l-a) 569, (s,+8S,) o,(s,) 


(1l-a) (1-a) 659, (s ‘e 


is erp 

aad (1-5) 6, (ys, +8S,) >, (8s)) op (s,) 
(1-a) ad (1-5) 6, (YS) 6, (8S}) >, (S,) 

a (1-a) 6 (1-5) $,. (ys, +885) >, (s,) 

(1-a) (1-a) 5 (1-6) 6, (ys,) o,(s,) 

a (1-5)5 9, (s,+8S,) >, (8S)) op (YSp) 
(1-a) a (1-6) 69. (S,) >, (88) o,(YS5) 
a(l-a) (1-6) do, (s,+8s,) >, (YS,) 

(1-a) (L-a) (1-5)5$, (81) >, (YS,) 

ao (1-5) (1-6) o, (ys, +8s5) 9, (88)) $2, (785) 
(1l-a)a(1l-s) (1-6) ¢,(ys)) $2 (8S,) op lysy) 


a(l~-a) (1-0) (1-95) $,(ys,+8S5) oR (YS5) 


(1-a) (1-a) (1-8) (1-8) 0, (78) og (785) 





(S)7S5) = [6o,(s,)+(1-6) o,(ys,)] 

x laason(s,+8S5)o,(8s)) + (1l-a)ado,(s,)o,(8s,) 

+ a(l-a)d¢,(s,+8s,) + (1l-a) (l-a)o¢,(s,) 

+ aa(l-o) d,(ys,t8s,) >, (8s)) 

+ (l-a)a (1-9) o,(7S8)) o,(8s,) + a(l-a) (1-5) 6, (ys) +8s,) 


7 (1-a) (1-a)(1-6) 6, (ys,)] 


(S)7S>) = [oon (Ss5)+(1-5) Oo, (ys,) | (Gite Cele) 


w) 
eae 1 


x ladon (8s) +(l-a) ]xfado,(s)+8s,)+(1l-a)5o,(s,) 


+ a(l-3)o,(ys,+8s,)+(1-a) (1-5) o,(¥s)) 1 


For the special cases of the EMA(1) process, III.C.7.1 reduces 


to the results given in Lawrence and Lewis [Ref. 5]. 


D. THE MOVING MINIMUM MODEL 
IL Gy ilinke ifereh bret sales el 
Another possible scheme that can be used to generate 
one-dependent sequences of random variables with marginal 
Exponential distribution is the so-called minimum model. With 
this model the {xX} Sequence is generated by taking the moving 
Minimum value of two Exponential random variables. The proposed 


generation scheme is 


An = MIN(E, ,DE_ feb. Deals 1) 


op 


oo 





where {Xt, me= 1,2,...}) is a sequence of random variables 
Peon marginal Exponential distribution, {Es ie Oe 1S 
an iid sequence of Exponential random variables with unit 
mean, and b > 0. This will produce an {X}} sequence with a 
b+1 b 

rate of 7 and an expected value of 5e1° This expected 
Mame produces one difficulty since it is a function of the 
parameter, b. This complicates comparisons between results 


with different parameter values and decreases the value of 


scatter plots and sample paths. However, this difficulty 


can be easily removed by multiplying the i by ol The 
generation scheme then becomes 
Sie b+] 
X = =p Xt = MIN(IA1E,,[b+1]E, 4), (Gieieis Dials) 


with ee and b defined as before. The {x} has a rate of 
one and, hence, an expected value of one. This facilitates 
comparisons for different parameter values with the NEMA(1) 
discussed in III.C which produces random variables with unit 
means. 

The investigation of the moving minimum model is moti- 
vated by the previous result in III.C.2 that linear additive 
models have a constrained range of serial correlation. The 
hope is that the non-linearity of the moving minimum model 
will obviate this constraint. The minimum scheme has been 
used by Tavares [Refs. 22 and 30] to generate first-order 
autoregressive exponential processes and by Marshall and 
Olkin [Ref. 31] to generate correlated bivariate Exponential 
variables. 


160 





Zeeeerrelation Structure 


The first-order serial correlation can be computed 


by the following approach 


ba! 


) = E((IMIN(I=—JE__,l 


D+1JE_) ] 


(MIN( ([2=*]E., [b+1]E,_,) 1). 


The terms inside the expected value can be made independent 


by conditioning on the value of E,: The E(x ) is then 


n¥17n 
Bouna by multiplying the conditional result by the density 


of En! and integrating. Implementing this approach we have 


es peal 
fe X) = pe eBay: (ty) 
b+1 - 
x(MIN([7E=ly, [b+1]E,_j1|E,_ Je Y ay 
<a ( (eturn( Ptlye [b+l]y) ]) Cli Drea) 
n+1*n 0 : b n+l’ ai i 
b+l 


x(B(MIN( [==ly, [b+1]E,_,) ])e 7ay 


The expected value of the minima can be calculated as follows: 


bx 
(b+l)y oss 
ap dl _ ee ©, Jsiqeul 
Ben (Eo, (b+1)y1) = J x (pay) e dx 
co aos 
+ ii (b+1) y (pey)e Ot lax 
(b+l)y 


Lisa 





-——-| (b+l)y 
b+1 b+1| 
Beni ——)e,,(b+l)yl) = -xe ; 
bx 
(b+l)y -—— = 
+ sa + (b+l)ye oy 
0 
bx 
(b+1l)y¥y ——— 
= =(b+l)ye by a my (— = b+la., 
Pee pen 
b+1 ee On. ea OY 
E(MIN[(—--)E_),(bt+l)y]) = es 2 ye Cit beDie 2) 
Samiiarly, 
b+] (b+l)y/b 4) “BRT 
B(MIN[(—=)y,(bt+1)E,_,]) = ; sles) dx 
oo eee 
USP Ie) 
| (bee (b+l)y/b -— 
_ b+1 Y b+1 
= -xe + 2 ax 
0 0 
5 pet) By 7 
(b+i9i7b eee 
eee aey. —-y/b ih, b+1 
_ a = ate (b+1) f (oa7) e dx 


0 


(bel y, _-y/. 
ae, 


G2 





joviwde = 1) (1 ewe Gare. D.2.3) 


E(MIN( (~==)y, (b+1)E,_41) 


ieeenG iiteb.2.2°and TIT.D.2.3 in III.D.2.1 produces 


b+l -y/b, - 
A, Xai) = j (PES) (1-e7 PY) (b+1) (1-e Y/?) e Yay 
D 2 
(bet) : see 2 jee — 
b(1+b+F) 
(b+1)* 
Sol 
Therefore, since E(x) =e 
oo b 
b“+b+l1 b°+b+1 
and 
CORR(X_ 4 /X,) = ee (egal eee ale) 
b°+b+1 
Thus the model allows a range of correlations from 
[0,51]. The minimum value is achieved when b is zero or in 


the limit as b tends to infinity. The maximum value is achieved 


wnen b is equal to one. An interesting aspect of the correla- 


Eron Structure of the moving minimum model is that reciprocal 


values of b produce equal correlations. This is a similar 


Kind of "invertibility" found for the other moving average 


models discussed in III.C. The range of b could be restricted 


IGS 





to the interval [0,1] without reducing the possible range of 
correlations. However, doing so, as with the NEMA(1) model, 
moulasrenonre the fact that in the non-normal case character- 
istics other than correlation may be different in the case of 
b and =. Also, as is the case for all the first-order moving 
average processes addressed in this paper, the correlation 
for lags greater than one are zero. So the study of correla- 
tion structure is limited to the study of the serial correla- 
tion with lag one. 
Eee cae tve eorrelation 

The range of possible correlations can be extended in 
a fashion similar to the NEMA(1) model (see III.C.2) by the 
use of correlated or antithetic variables. Using this approach 
the generation formula becomes for antithetic variables 


Orr ib 


eee one) (pt1) e! | (Siete eat) 


ys 


Peal variables are defined as in III.D.1.1 and {Ei, 


m= 0,1,...} is generated from the Eee sequence using the 
an 


» relationship EA = -1n (l-e Note that this implies that 


{Ets is also iid Exponential with unit mean. Again 
joe aa 


) = E[(MIN[(—-) E_, (b+1)E!_4]) 


loeaak 


x (MIN [ (—=—)E (b+1)E.] )] 


nei 
ema, CONdie1Oning on the value of EA! multiplying by the density 
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of En! enemwluntegrating produces 


oO 


zi ee 
ek) a pee inp lyr ier iE). ) 1 
b+1 ‘ al vod 
[MIN((---]E Lj. Jerid) IB = ic dy 
m(X xX ) ( ceimin( 2d [b+1]E' ,)]) CIE spas) 
n n+l bY! n-1 re 
b+l 


(E [MIN ( [===] E_ ,+~{b+L]1n[1-e"%]) ])e %ay 


aa 
Mae first expected value is identical to III.D.2.3. Thus 


E(MIN((224+)y,(b+1)E"_,]) = (b+1) (1-e"¥/®) Omri) 8 3) 


The second can be calculated as before. 


b+1, . _ We hs 
E(MIN{ (==) £, |, ,-(b+1)in(l-e *)]) 
-(b+1)1ln(l-e *) oes 
= ae ae 
baer 
0 
xo _ Dx 
6 i -(b+1)1ln(l-e “je te 
=(ped)intl—-e ~ ) 
~ Den. pey. 7 ey nox 
pet|-(btl)in(l-e 7) 4, (b+1)In(l-e"*) , 1 
“= b? ‘b+I 
0 0 
: ey 
eee in(ioe ye ne) 
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E(MIN( (P=) E_ |, ,-(b+1) In(1-e7¥) J) 


wey Le 

= (btl)In(l-e Ye? Me ©) 4 BEL) (y--Pinfl-e “1, 
= mee cy 
iorsmin(t=er! jer ° 


orp il 


E(MIN[(===)E Ode Slat (ler smeay ip 


need! oy 


ep b 
b 


b 


eee g = [j=e 2 ]-) (III.D.3.4) 


Peet emerng If Pes .3 and ITi.D.3.4 into III.D.3.2 yields 


7 tol -y/b, ,b+l - - 
E(X Xj) = f (Pet) (re ¥/P) (BES) (1-[1-e7¥ 1") e Yay 
b+1 : 
= (*E)*s eYay- AEs Aye % “ay 
0 0 
(PEA) Ay oY (1-27) Pay 
0 
poe 2 -(14E)y -y.b 
=F aw) J = (l-e Y) dy 
0 


The first two integrals are trivial. In the third the change 


of variable z= Cice dz = e ~dy makes that integral straight- 


forward. In the last integral, the change of variable u = e?, 
8 = dy makes that integral recognizable as the integral of 


a Beta random variable. Using these changes of variables and 
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making the appropriate changes to the limits of integration 


produces 


i 0 
b+1,2 a. Ga Pee 


Be tL! ee Sale "pie ) b 


b+l b+l b+l i fe 


= 7S) - Ce ens are) SKE 


It (145) lia) 


 (24b+e) 
enter & 
pt (6) br (b) 
1 1 i 
(1l+b+-) (b+=) Ip (b+=) 


2 
bol (s) TP (b) 


= oS a a, es 
(b +b+1) (b°+1) I (b+) 


Then 


yee | 
BT (s) 1 (b) 


(b*+b+1) (b°+1) I (b+) 
and 


brisjr tb) 
_ Se a Ss | (ITI.D.3.5) 


eT? 5 5 


CORR(X, ,X 7 
(b“+b+1) (b°+1) I (b+=) 


Like the expression for positive correlation, this 
expression is also symmetric with respect to reciprocal values 
of the parameter. It attains a minimum value of minus one- 


third when the parameter value is one. Graphs of the 


IL 





correlation as a function of the parameter b for both posi- 

tive and negative correlations are provided in Figures III.D.3.1 
and III.D.3.2, respectively. Unlike the NEMA(1) model which 
requires an additional parameter (see equation III.B.6) to 
achieve a full range of negative correlations, the moving 
minimum model can achieve its full range with the single 
parameter b and an antithetic sequence {Et}. 


wemwoint Density of x and x 


+] 
Calculation of the joint density of x and Xe] is 
possible using a conditioning argument to determine 
eee =| 2 = 2) and 9 (aS Zz). These values along with 


emewprObdoillity that EW takes on a given range of values are 
sufficient to determine the joint distribution function of XS 
and X44: The form of the distribution will vary depending on 
whether one is above or below the line X.,, = bxX,. The joint 
density, where it exists, is determined by differentiating 

Pie aistribution function. 


From III.D.1.2 we nave 


_ ete] 
sa = MIN ([-=>-]E,, [b+1]E, 1) 
Then 
( l if (2f*)z < x, 
P(X. <x|/E_ =z) = CT teeta) 
n=— ig) a p 
| tos 27k ae (2B) z 5 oe. 
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The first result in the above is obvious. To justify the 


second, consider 


= Dx = 
eee = 2 > 7) = P({b+1lJE | <x) 
a < 
ine eee pasa 
eee 
a b+1 
; _ b+l 
Since | = MIN(([=>-JEL,,, (bt+llE,), then 
oo 
( 1 Ne eo Srey 
eee SY |B = 2) : _» 
b+l . y 
l-e ie. "2 > DiL° 


Remeensiader the Joint distribution, note 
xb 4 = 
ma < 5+ (i.e. when you are above the line bx = 
of possible z values can be broken up into three 


mecgure IilT.D.4.1. Then, since ES 1s Exponential 


mean, 


m = bx 

P({z <¢« REGION i) = Piz rae ee 
__ Dx 

Piz e REGION 1) = 1 aoe 


Piz <¢ REGION 2) 


" 


OT ip 34a) 


that when 
y), the range 
regions. See 


Wacacie) (bhenise 


Gr ae Oy by re, 





ee Tails) ie 


L+¢ L+4 
en ee eR hs B eN 


E NOIOSYd c NOITOSH | NUILOsd 








PS ese ae etd ie 


| + L+4 
of pe x4 cme 


ee .o° -o —~- .  eeee e 


CF NOIOIH' 2 NOISFH' T NOLOGY 
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EZ eR EGLON 2) 


P(z € REGION 3) 


Piz «© REGION 3) 


Now bv definition 


iI 


il 


b+l _ 
oe 
NZ er! 
2 a 
b+lL 


ieee kee Sy|2 ¢« REGION i) = 


Pix 
n 


because when conditioned on the value of En! these probabili- 


ties are independent. 


mee b.4,2, and the definition of the regions in Figure III.D.4. 


P(X <x,X,1<ylz 


Ee en SY |Z 


P(X <x,X,,, <y|z ¢ REGION 3) 


Prniaeeremnecmlts Of P1i-D.4.3 through IIL.D.4.5 and III.D.4./7 


SauouGh Lhinnna.o we Can Comoute the joint distribution of Xx, 


and SL 


~ 


= 


— 


+1 


REGION 1) 


REGION 2) 


i 


P(X <xlz 
n— 


<y|z ¢« REGION i) 


when y > bx by uSing the relation 


(Tia. Dees 4) 


(LEB) 


¢€ REGION i) 


Using the above equation, 


(i ieee 4.. 


te neab< 42h; 


ia ee Dye 


(TRE eb. 


6 ) 


ao) 


0) 


9) 








Il 4 Ww 
td 
Ps 


P(X <x,X1) <y) <x,X,,, <y|2 « REGION i) fel ate IkO)) 


x P(zZ « REGION 1) 


Dx x Dx 
a7 gee b+ly se b+1) | Shae 18 a ay 
_ xX _ by 
saa (oS 2 alee oe Ve ENO 
DJs Dx y nea. ey 
= je Ptl Seri Daa. eee oo b+l | B+] 
meee =" 
b+l -y (Sar) (xt [btlly) 
co — ea a = 
ry 
ae eee loreal 
P(X, <x,X.,,<y) = l-e “-e 7+e (ii Teer ay veel) 


Similarly, when y < bx (i.e. when you are below the 
line bx = y), the range of possible z values can be separated 


miqco three regions. See Figure III.D.4.2. Then 


P(z e€ REGION 1) = Pe Set 
a oe 
b+L 
Piz © REGION 2) = l1-e ; (iid sD -2ese) 
2 4 bx 
P(Z « REGION 2) Play < 25 p47)" 
> a __bx_ 
P(z ¢ REGION 2) = e PTL | Pl. (Taba aiken 


EA 


aa 


REGION 3) = P(z > 





P(z « REGION 3) 


Pace mir. D.4.), IZTI.p.4.2, IfI.p.4. 
ene regions in Figure III.D.4.2, 


for the given region. 


P(X, <X,X,,, <y|2 ¢ REGION 1) = 
P(X <x,X_,, <yiZ « REGION 2) = 
pe oy Kea Sy 2 - REGION 3) = 
~_*_ 
Ka(iaet es ) 


Pemolnang LILT.D.4.10 with III.D.4.12 


(Ciaier, Did lie.) 


6, and the definitions of 


the following results hold 


i Gar pe ons) 
_by 

eee (ici pees 
ae 

(ise 25>) (TTT. pet 7) 


EaEouUgM “LE D.s. 17 yLrelds 


wer bx < y 
eo ae 
P(X, <%,X,,,<y) = l-e Ye *te 97? (III.D.4.18) 
x 
Let Ft (x,y) = (ee ee Ememaistrlbue1Ton 
Bunction of Xo and K+] when bx < y; and let + (x,y) be the 
joint density of x and K+] when bx < y. Then 
—T-¥ 
finy) = PE Pixy) = 3 fire te Pte OI 
_—_ 
= te Ye Joka * 


Lvs 





f(x,y) = (gpyle bx < y; x > 0. Goes ale) 


For Dx 1< Vener atstrt button function, P“ (x,y), is 
=) _ 3 oa 2 


l-e *7-e “+e melt oe years the jOolnt density of se 


and X +] when bx > y, then 


2 oh, 2 ae -y -x -x- pty 
fe(x,y) = 3x Dy F°(x,y) = ae tee -e “+e ] 
ee) i 
3 - b b+1 
= <-fe Y-() 
aX b+1 
2 b -x-p dy 
1s ay) = (57) Sy DIO LGh yY > 0 


Note that there is a positive probability that the 
point (Xe X47) lies on the line bx = y. This probability 


can be computed as follows. We have 





- b+1l 
x = MIN([—>-]E,, [b+1L]E__) 
x b+l1 
x4] = MIN ( laa=eee Saree [b+1]E) 
mme point (X ,xX ) lies on the line bx = when X = (tl). 
P jae ie acl J se p 
and X +1 = (b+1)E. Now 
_ -b+l1 : _ 
P(X = Sp I ERtXy4, = Pet lE,) 
_ jase |: ; lose It 
= P([--1B, < [btlJE 7 (etllE, < (=[-I1E,44) 
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The events in the right hand side can be made independent by 


Concitiomang on the value of En! Then 


ey Spl ; e 
P(X = (=—1E? X45] = [btlle,) 


a losp il b+1 E -y 
ay ee he elo eee (>ellye< [|e |E_ =v)e “dy 


Lo @) 


y 
a1 ieee 


1 > by)e *dy 





+> 
1 = 1 -(btl+a)y 
= I ii (b+1+F) = dy 
b+l+— 0 
b 
b+1 b 
P(X_ = [=—] E_;X =[bt+1l]E_) = —~—. Os eee le) 
n b n’~ n+l n b-+b4] 


Because there is a positive probability of lying on the line 
bx = y, the moving minimum model can be said to have a line 
degeneracy. An important implication of the positive proba- 
Bality of (XX 47) lying on the line bx = y is that the 
moving minimum model will produce runs of values of constant 
ratio b. The values of ey in these runs will be decreasing, 
equal, or increasing for b less than, equal to, or greater 


than one, respectively. The length of the runs will be geo- 


metrically distributed with parameter —2— for the positive 
bet b- L 
correlation case. It was this kind of degeneracy in the 


Exponential autoregressive model, EAR(1), that proved to be 


one of the model's major weaknesses. The degeneracy also 
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occurs in the Tavares autoregressive model and the bivariate 
Exponential pairs derived by Marshall and Olkin. 

The probability of lying above or below the line bx = y 
can be easily found by integrating the appropriate joint den- 


peyeover the area desired. Thus, for bx < y, 


P(lying above bx=y) f f £* (x,y) dydx 
0 bx 


J f (ye 
) ee OS 


P(lying above bx=y) = -—-==--——— (Libap 4222) 


Semvlarly for bx > y 


oO 


jop:< 
P(lying below beag= yy) =" f £* (x,y) dydx 
0 0 


| 
—, 
| 
% 
7 
i 
(9) 
oO 
4- 
oe 
Qu 
~ 
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PeeeeenoaitiOonal Expectation and P(X Sen) 
i eee a “+ | n 


Besides the correlation coefficient, there are two 
Sener. Characterizations of the joint distribution of x and 


K+] which we have considered. They are the conditional 


expectations and the P(X, ORE) Both of these quantities 


+] 


can be derived by considering the four possible sets of values 


Or x and X computing the probability of each set occurring, 


piace 
and weighting the conditional expectation or probability by 
its probability of occurrence. 


Pir stt he sorObability of occurrence for each set of 


values must be calculated. Consider the case where XA = OFS) E. 
b+1 
eee 85 Rpt 
b+l1 — =.oReth 
ee En’ Sn = | b JE 42? 


b+l b+l1 


= a [b+1]E i, [I= ae [b+1 JE) 


hic 
By conditioning on the value of Eye the events on the RHS 
become independent. The calculation then proceeds in a 
moetraightforward way. 
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b+1 - 
The second case when X = (=—)E_ and X = (bt+tl)E has 
¢l b n n+l i 


already been computed. The result is III.D.4.21 and is re- 
peated here as 
b+1 


en fp Ea Xn 


ie [b+1]E) = eee (aiviebe D5.) 


The third case is when X, = (btl)E__ 


Here we proceed as in the first case. 
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The final case is when x = (b+1) Ey and | = (b+1)E.. As 


before 
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nN 


‘a b+1 b+1 yar 
Set 1E 23 < ee ogi Mee EO ee yiie ~ dy 
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The conditional expectations can now be written by 


inspection. 
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(PE Soar = FEI = 
a (PRRIE Xi) = [b+lIE, PY 

en - (DALIE, 1X4. = CEE — 
pee- (beU)E! 4, X,,, = (btlle._, b+1 


Weighting these conditional expectations by the probabilities 


Meelis oases! through LLL.D.5.4 yields the final result. 


4 
B(x, .,/%,=yY) = pei 4), ~ Bree 1)P(case i) 
b+ ae b 
= (2) (2, 1 + (by) (5 —) 
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Seal 
B(X_,{X,=y) = 1+ —+=. (III.D.5.5) 
ee b“+b+l 


t is quite surprising that the regression of Xney On x is 
linear in y, considering the non-linearity of the process which 
generates the er - 

The conditional expectation of x given Been can be 


derived with equal dispatch. 


Eo 2 


_ 


- 





<aSe B(X) 1X4) =y) 


<n a ae ope) a 
X= PEIE X41 = (btlle = 
B= (PHB, 1X4) = PEE a 
X= IbtllJE, _ 1X4) = [btliz, b+1 


Bomeg LEE. b.5.1 through III.D.5.4 as before 


4 
mx |X.) =y) = ae atv case i)P (case i) 
b+l b> y b 
i. ia) + ) (—>—) 
[b+1] [b°+b+1] b° +b+l 
“ (b+L) (x2) + (b+1) (=, —_) 
b~+b+l [b+1] [b-+b+1] 
| _btty 
Ex | xX = y) = + 1. (iD. 5260) 
i +1 b<+b+41 


The probability that mal 1s greater than x, can also 
be easily computed if one is careful to differentiate between 


mpmeecase where 6 < 1 and b> il. 


Case P(X412-4%,) 


NO) fl oe 


os 





Gmeeete Ds 1, 
om z PIE XY ~ [b+1JEL | 5 
l1 if bol 
“4 se Sarl il 
R= Wet En)  Xa4y = EH Ena bi 
Me pti je. .,X .. = (b+1]E * 
n n-1’ “n+l n 2 


iipcewhenwia < | we have, using III.D.5.1 through III.D.5.4 


4 
fom “A = pee aaa * Xl ease i) P(case i) 
3 
Fy) 
(iscpll || si arene il | b +b+l 
+ () (———+ 
(b+1] [b +b+1] 
JP) OR ax) = iL =. i ee itp ape (Lik Deon) 
ee 2 “ (bt) (674541) = 
Meerniimlan COMpucatlon with b > 1 again using III.D.5.1 
Mmenrough III.D.5.4 yields 
So ae b 
eee x) = 54 ja 2 dL Gault D aoc) 


Goel (GC -use a 


Thus a graph of Eee aa) will have a discontinuity at b= 1 
when case 2 switches from a probability of zero to a probability 


of one. This graph is presented as Figure III.D.5.1. The 
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minimum value of one-third occurs at b = 1. The maximum value 
of two-thirds occurs at b = i The moving minimum model, 


therefore, has a greater range of values for the Pit 7 Xp? 
than does the NEMA(1) model. However, the greater range for 
the B(x +1 7 Xp) and greater range of correlations must be 
balanced against the degeneracy of the model. 

As was noted with the NEMA(1) model, the correlation 
in non-normal models does not define the joint properties of 
ee and Xnel° Although the cases of b and = are indistinguish- 
able from the viewpoint of correlation (see III.D.2.4 and 
III.D.3.5), these cases will have significantly different 
sample paths as indicated by III.D.5.7, III.D.5.8, and the 
meseusslon Of runs up and down in III.D.4. 

Three examples of sample paths for different b values 
are given in Figures III.D.5.2 through III.D.5.4. The degen- 
eracy of the model is clearly present in the sample paths as 
a tendency to produce runs of equal, increasing or decreasing 
values, respectively. A comparison of Figures III.D.5.3 and 
ITII.D.5.4 quickly demonstrates that while these two sample 
- paths have the same correlation, they produce significantly 
different Osa! Sequences. This is a graphic indication that 
non-normal processes are not determined solely by their 
Ssorrelativon structure. 

Paguees DTTI1-D.5.5 through ILI1.D.5./7 are the scatter 


plots associated with the sample paths in Figures III.D.5.2 


through III.D.5.4, respectively. Here, too, the degeneracy 
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of the model is clearly present in the tendency of the 


(X 7X plots to lie on the line XO = bX. The slope 


+] 


of this line determines whether the runs are of equal, in- 


sual 


creasing or decreasing value. 


6. Conditional Expectation and P(X 412% Wr ror 
Antithetic Variables se 


Results similar to those obtained in III.D.5 can be 
obtained for the moving minimum model with negative correlation. 
The procedure for determining the conditional expectations and 


the probability that x, 1s greater than Xo MSang -antLthecic 


+] 
variables is exactly the same as that in the previous section. 
First, the probability of each of the four possible combina- 


Erons of Xo and xX values is computed, the conditional 


+] 
expectation or probability is computed for each case, and 

the final weighted sum of conditional expectations or proba- 
bilities is finally computed. In one instance no closed form 
answer is available and numerical procedures are used. 


Recall that the generation scheme when using antithetic 


Variables is 


= b+1 ' 
x = MIN([P=JE_,(b+1]E!_j], (III.D.6.1) 
where {E n= 0,l,...} is an iid sequence of Exponentially 
distributed random variables with unit mean, {Eos n= 0, lee 


is generated from the {EL} sequence by the relationship 
15) 
Ee = — In {1l-e ") which implies that {E,} is also iid Exponen- 


Piabewhen Unit Mean, b > 0. 


ToS 





bl 


First, consider the case where Xo = ae oe and 
_ -bt+l 
a ST Ene: Then 
_ -b+l1 _ ,otl 
a, = | b ene b 8 na 
b+1 ' b+1 


= ee USreeee eoinens [b+1L]JE 1’ ee Ea = [b+1]E*) 


n= 
Using the standard conditioning argument produces 


_ ponedk > 
mo = I b ce a = | b eae 


Datel 
b 


JE] <7 tb+llin(l-e 7) JE. =y)e “ay 


f PUPEAly < tb+11E,, | a 


jp (EY, > H) P(E.) <~bin[l-e *1)e ay 


f wee ae fier |) e ay 


0 
1 uN 
eee) 0 = _ tiiereS} 
= fe P ay - f{ (ne %)Pte%) P ay 
0 


The first integral is straightforward. In the second integral, 
the change of variable u = e?, <3 = dy makes this integral 
recognizable as the integral of a Beta random variable. Making 


this change of variable and making the appropriate changes in 


m@e limes Of integration produces 
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+1) 


Oo |r 


du 


L+b)P (145) 


ut 
[ (2+b+F) 


Ae) [b+1]E. 


n ai! 


b°T (b) T (2) 


*+b+1) (b*+1)P (b+) 


ehgk 


(EEL. 2D ao) 


In the second caSe, x = IIE, and K+] = [b+1]E). Pro- 


ceeding as before 


Alea 
P(X, = | b JEW X, 


lorie 


P ( ae! EA Ss [b+1]E)_y, [b+1] EF < ae = 


oO 


f PC 
0 


oe) 


local 


b 


a [b+ 1]E) 


ipa 


falar lL 
este EO Ss i>eljinit—e9-] < | 


Vv 


' a - ee = 
J PCE a1 >) P (Ene? bin[l-e *])e “dy 


f ee een en dy 
0 


~ (145) y 


{ (-e7%)Pe ay 
0 


lithe Be 


ogere 
b 





= ay 
Pee a ay 





This is the same integral as the second integral in the first 


case. Thus 


bel 


PAX, = (pm E ys X47 = [etl IED) 
b°r (b) F (=) 
= (Tit sp oes) 
omit) (b +1) T (b+5) 

Nexc Consider the case where X= (bt+1)E!_, and 

aor 
on = (“) Ee Then 

b+1 


eee re pen ~ ETE 


. be orl 
BK [b+1]E* lee, laa B 


4) < tbt1lE! 


cC 


jf P((b+LI Ey < (Fly, PE JE, ,, < (btl]in{l-e ¥] |g, =y)e Yay 





= f aS ee) a -bln{l-e %])e *dy 
0 
& 3L 
oo fee) — lt) 7 co : oe : - aiee=) 
= feYay-fe  ey-f e¥(-e ¥)Pay+f (r-e%)%e%) ay 


The first two integrals do not present a problem. Making the 
change of variable z = jane az = e Ydy makes the third inte- 
gral easy. The last integral is the same as the second integral 


in the first case. So 
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bai 


ea’ “41 = | b IE ny? 
: ‘ po (b) P(e) 
sels pa * ae 
a (b“+b+1) (b“+1) I (b+=) 
| ese 
P(X, Ee Sn. | 5 Baye? 


b*r (b) r (2) 


(b?+b+1) (b741)T (b+2) 


Finally, consider the case when Xx = (b+1) ET 
—_ t 
X41 = (D+1)E.. Then 
= t = t 
P(X, [b+1) EO Xana [b+1]E*) 
b+l1 b+1l 


P ( [b+1]E' < i En! [b+1]JEX < Fed eo 


lo) 


{ P((b+llet< PE 
0 


= ly,-{b+l]lin[{l-e 7] < [ 





oO 


i Ve = em a 
j PCE n 1 <6) Plena bln{l-e *])e 7dy 


Lo, @) 


{ (1-e7*) 
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b- 
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Cw lem ete.) 
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(Tit. pe6 a4) 


and 
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bel 
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JE ad |E_ =y) e Yay 


These integrals are the same as the third and fourth integrals 


in case three. Therefore, 


doe 





P(X, = [b+1]E'_,,X_,, = (b+1]E") 


sgl 


1 bP (b) r (2) 


= 1 aes aie oe (LiT.D.645) 


Pt" (b*+b+1) (b241)T (b+) 
The conditional expectations given a specific case 


of the values of x and x can be written by inspection. 


+1 
Hence, 
ease MOR ee 

_ -btl _ .btl b+l 

ee epee nl nl | | b JE nad! b 
Lee 

_ woerdl _ ” - oS pk 
= [ 5 JE Xa) = [b+1 JE; (o+1)ln(l-e ie 

- . _ bel. b+l 
se a Xn bp Ens? D 
X= [b+llE!_)/X, ) = (btllet: b+l. 


Combining these results with equations III.D.6.2 through 
fier. b.6.5 and letting 


bol (b) T (=) 


(b*+b+1) (b7+1) f (b+E) 


we have 
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- 


Bex |X, = y) = (2 ee Pease 1)P(case i) (2 Ey. 6280) 
PY 
= (2-6) ®E) - (pti) in(i-e Pt?) 
b+1 1 
+ G(PEA) + (93-6) (b+) 
py. 
Peo, |x =y) = 2-G(btl) (1+In{1-e °°") cep ce 


Similarly, we can derive the expression for E(x |X.) =y) 


Case BAe Xo =) 
b+1 b+l1 b+l1 
oe b JE X47 = | b JE nei bi 
b+1 b+1 “b+tl 
= (arom EA +1 = [b+1 JE; a a in(l-e ie 
a _ bth 
Rebel IEY | Xo.) = (NEL? Coe / 
= ' = Lee 
Xo [b+] EB! Xa] [b+1]E}; earls 


iaeneising £1L.D.6.2 through III.D.6.5 and again letting 





b*r (b) P(e) 
= = 2 i 
(b~+b+1) (b e41)r (b+) 
4 
E(X |X.) =y) = ee a nei renee 1)P(case 1) 
at Rea > 
= (—— ) + (b+) Gt (b+1) (S45 G) 





noo 





b+l ~B+1 
BO | Xa = Y) 2- eC .) (1 + in[l-e 


The probability that x is greater than x can be 


+1 


approached in the same fashion as the conditional expecta- 


tions. The second case will be reserved for individual 


attention. 
Case P(X +1>X_) 
Sees ——n———s n> 
een we _ ,-btl it 
te pr ne ned op Ene 2 
- _ -bt+l 1 
X= [btllE, 1X41 = b JE nad b+l 
X = [b+1]E' ,,X_,, = [b+1]E! i 
n n-1’ “n+l n 2 
The second case, X = (Pte and X = [b+lL]E', does not 
i. 1 b n n+l Ee 


allow a closed form solution. We get 


Dtlip 


eee) = P({b+1]E, > ae e 


E Bee 
ie ie — ae) 


= P(-[b+l]1ln[l-e 5 _ 


> [ 


iG 
= © lise eee mn) > =>) 
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Px 


“YO. b “Yo 


Menee, if we can find Yo such that (l-e )~ =e Pree tl 


y 
P(X 0 


a = l-e The required solution can be found by 


n+l 
numerical means for any given value of b. A computer program 


6 


to determine y, to an accuracy of MO forvagiven value of 


b and to compute P(X, = OS was prepared. A graph of the 


ale 
results is presented as Figure III.D.6.1. When using anti- 
thetic variables, the moving minimum model has a restricted 


>X _). The maximum value 


range of possible values for P(X41 bs 


of approximately 0.509 occurs at about 0.30. The minimum 
value of approximately 0.491 occurs at about 3.33. 

This small range of values for the Exe is 
shown in the relative indistinguishability among the sample 
PeeoeGisolayea in Figures IIL.D.6.2 through III.D.6.4. Of 
more interest are the scatter plots presented in Figures III.D.6.5 
through III.D.6.7. In these plots the degeneracy of the moving 
Minimum model is clearly displayed. When xe achieves a value 


of y based on Ea! then Xo is constrained to have a value less 


= 
than -ln(l-e ”). 


+] 


In the case where equality is achieved, the 


second case in the discussion in this section, the point (Xo, 


—% 2 sc] 


) lies on the curve e ee —MiPeeiits constraint 


n+l 
is clearly evident in the scatter plots. Thus the moving 
Minimum model displays a degenerate behavior for negative 


correlations as well as positive. 
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E. THE BETA-EXPONENTIAL MODEL 
fe ee Oduct lon 
A third method that can be used to generate correlated, 
Marginally Exponentially distributed random variables is a 
special case of the Beta-Gamma model given in Lawrance and 
Lewis [Ref. 18]. This model generates an {x} sequence using 


mac relation 


xX, = 8, (a-1l-qEL + Bl -d GEL _y Bvt lic: Coren ea 7 ee yy 
where (B.(q,l-q), ie—-  otomaneiid sequence ios Beta 

random variables, {B_(1-q,q), n= 21;72,...! is an iid sequence 

of Beta random variables, {Eve Ce) sma scouenec 


of Exponential random variables with unit mean, {B_(q,l-q)}, 
{B_(1-q,q)}, and a are mutually independent, and 0 < q <1. 


The density for a Beta(m,n) variable is 


Tr (m+n) m— 1 


aaa , : 
Pon) PGE X (1-x) @) i x Ss i ne 0; n> Oe Git B12) 


.iIn practice the Beta random variables can be generated from 


two Gamma distributed random variables using the relationship 


G (m) 
G(m) +G (n) 


1s a Gamma random variable with a shape parameter of K and 


B(m,n) = from Kotz and Johnson [Ref. 19], where G(K) 
a scale parameter of one. 

This is a special case of the Gamma model considered 
in Chapter II of this thesis. It works because,as described 
Syebewls [Ref. L0],in IIL.E.1.1 the product of the B_(q,1-q) 
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variable and the EL variable is a Gamma (q) variable. Simi- 
ila the product of the BL (4-4,4q) and the E-l variables is 
a Gamma (l-q), independent of the Gamma (q) variable. Conse- 
quently, their sum is an Exponential variable, X,: The {X} 
process is clearly one-dependent, as for the NEMA(1) process. 

Because of a lack of a closed form solution for the 
integral of the Beta density when the limits of integration 
are not from zero to one, this model is the least tractable 
of those considered in Chapter III of this thesis. However, 
its correlation structure can be determined, an expression 
for the Laplace-Stieltjes transform of a sum of r random varia- 
bles can be derived, and the probability of X41 being greater 
than Ke can be simulated. An advantage of this model is that 
1t extends directly to moving average Gamma processes (see 
@iapter Lf). This extension is not possible with the NEMA(1) 
or the moving minimum model. 

2. Correlation Structure, Positive and Negative 
The correlation structure can be determined uSing a 


standard approach. We have using III.E.1.1 


XiXna, = (By (a ta) ELtB, (1-q,q) EL_)) (8,4) (ge 1-Q EL 
i yay ah Same qe, 
~ arene eee) Rect ed) Boy = ne Peas Set eam GONE ae el eee bay eal 


4 B41 (i> G; q) BL (Gal= -q) E+ +B ne] (1794) BL (1-d,q) BLE _y 
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Taking expectations and using the iid nature and independence 


ef the {B_(q,1-q) }, {B_(1-q,q)}, and saa! sequences yields 


2 
E(X,X.4)) = q* + a(l-q) + 2q(1-a) + (1-q)* 
Hence, 
COV(X,, X14) =e ci) 
and 
CORR(X,,X_,4) =eeql—e) 5 0 '< gq <.t. (Cite) 


As with the other linear additive models, this correlation is 
double valued and lies in the range (0,2). 

The range of possible correlations can be extended 
to negative values by modifying the generation formula by 


including an {Et} sequence. Thus 


= = eel t 
XS Boi(qr1 qQ)E, + BL (1l-q,q)E,_1, (Ai <2 a2) 
where all random variables and constants are as defined for 
Mmeber. tl. and Br n= 0,l1,...} is an iid sequence having a 
specified correlation with the ee seguence. In particular, 

En and ES may be an antithetic pair. The correlation of the 


SS) using II.E.2.2 can be determined in the same way as before. 


Consequently, 


Pan al 





Ke] = (BL (q,1-q)E, +B (1-q,q)E!_)) 


(B+) $4 1-g) Fy + By, (1-4 q) ED) 
mete de) ce ea) Be Boag $a 1-g) BL (-a,d) 24 Bea 


+ BL 4, (1-4/9) B, (q,1-q) ELE? +B.) (1-a,q) BL (1-q,q) EXE? _| 


Taking expectations as before 


I) 2 
E(X X41) = @ tall-q) +q(1-q) [COV(E,,E*) +1]}+(1-q) 
= — t 
di eae) ( 1 gq) COV(E,,E,) 
Therefore, 
= ’ 
COV(X_,X_,4) qu q) COV(E_,E_) 
and 
: _ 7 
CORR(X_,X_ 44) = gq(l q)CORR(E_,E!), OD < cps is Cae ER ee) 
When ES = EA! Phe veemaalatien ts Onesand IIL.E.2.3 


meaduces £O ITT.E.2.2. When En and Ee are antithetic the 
correlation is -0.6449 and negative correlations result. When 
gq is 0.50, the correlation for the {X,} sequence falls in the 
(-0.16,0.25) range depending on the correlation between E 


ene &.". 
n 
MN 








a Haplace=stieltijes Transform of a Sum 


sg 
vk — ) > Ge Os 2 ree a 


where {X, } agmemactined by PiTsE.1.1. Then 


+ ... + B(q,1-q)E +B (1-q,q)E,_ 


B(q,1-q) EB +B, (1-q,q) E)+[B, (q,1-q) +B, (1-q,q) JE, 


+ ...{[B_, (q,1-q) +B, (1-q,q)JE__, 


eeyat 


Let dm = E(e ‘). Then using the iid nature and independence 
a 


of {B.(q,1-q)}, {B_(1-q,q}, and {E.} 


-s(B_(q,1-q)E_+B, (1-q,q) E)+{B, (q,1-q) +B, (1-q,q) IE, 


+...+{B__, (q,1-q)+B_(1-q,q) }E__,] 
isa = ile ial 6 r-l 


-SsB_(q,1l-q)E_ ~sB, (1-q,q)E, 


= E(e )E(e ) 


-s[(B, (q,1l-q) +B, (1-q,q]E, 
x E(e ) 


cee —cyra. (l—q,q))E 
ae Og = r~i @ r-1 


ONS 





op (8) = Uppgl Teg egl yg ts) 1 
tre! iy, (TO 
where 
a alls et eee! |* 5, 


The Laplace transform of the sum of two Beta random variables 
is a confluent hypergeometric function. Its form is too compli- 
cated to be of significant value in deriving the analytic 
behavior of the Beta-Exponential model. 
4. Empirical SOS pearsay 
Because of the presence of the Beta random variables, 


eae Probability of x being greater can not be analytically 


sal 
determined with a reasonable amount of effort. In an attempt 
to establish a range for this probability, a simulation was 


used. In order to achieve a precision of at least 0.001, 


Sixty-eight thousand comparisons were generated for each of 


- ten random number seeds. The Beta random variables were 


generated using the Kotz and Johnson [{Ref. 32] relation 


= G(m) ; B | 
B(m,n) = eimeeiay S<plarnec im fil Esl. The Exponential 


sequences were generated by a call to a standard generator of 
Exponentials. When the simulation was run for nineteen values 


femegecrom 0.05 to 0.95 in steps of 0.05, the P(X, ak) was 


+1] 
Heev0 for all values of q. 
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+] cael 


greater than xX is constant at a value of 0.500, reminiscent 


Although the empirical probability that Xo 


of the autoregressive model of Chapter II.B.6, the distribu- 
moon Of RE oes 1s not symmetric and no simple proof for 
this situation has been found. 

The low serial correlation and the apparent invaria- 
Serity of the ee makes the use of sample paths and 
scatter plots of little value. Samples are provided in Figures 


meine. t Enrough LTIL.E.4.1l2. 


Ze 





“OOT 


to A tit sano rs 


dg SEWN 
00°Sé 0S°29 00°0S 


I i 


XAQN | 


6h’ O? CIN) X<4T+N) Xi d 
81°O *OHY JTdWHsS 
61°C #OHY FNL 
S2°0:0 
HLlod 4a ldWbS blad 


OSiadeee OO Sie BS ce 00° 


Te 


ZS 





OCT 


ab 


US ..3 


Go Gls NAb elslOhoa hcl 


HAIGWAN XJAOUNI 


a 

BHO C(N) X< (14N) X) d os 
he°O *OHY J1dWbS 

22°0 *OHY JNHL : 

19°O?0 ; 


HLGbd aldWus tblad 


Ut 


ONS | 





“OOT 


OS 28 


Sige Gl tile Gholi 0a) ins! 


HAGWAN XAONI 


GUE Se OSme9 0005 BIS) (Ete my Th 


1°? ((N) X< (T4+N) XJ) ¢ 
he°O #OHY J TdWYS 
ec°O #OHYH JNU 
ce°or0 
HLoGd 4a ldWbS blad 


OO'E 
J lduldua 


OS‘h 
Slamien 


Z2arS 





"BE TF SCHT TER -PLOT 
W:Q.25 
Peewee fs Onl dg 


moet. Ano: 0.18 
O ‘EXTREME VALUES NOT PLOTTED 


6.40 





FIGURE III.E.4.4 


Deeg 





“BETA SCATTER PLOT 


We ex) 


Gel, ey 
TRU taal = Cle 2 


— SAMPLE RHO: O.2u 
O EXTREME WG le S = NG) ele Gee 


.? 

4 

my * 

6 

* ve 

e IL. e 
e J é O% 
nye s oie o. 
© se 3° ” 
co! 

ef a e,e 
cin e ae 
ee: A 
eee. ae iS 
e een? 

. ry no 





ERGO Lite .4..5 


Z20 





UO 


BO 


"BETA SCATTER PLOT 
BS 6.0133 
TRUE RHO: 0.22 


- SAMPLE RHO: OQ. 24 
3 EXTREME VALUES NOT PLOTTED 





EIGURE G2 lise ad 26 


PAM 





“OOT 


Loe Aster dane ta 


d4IGWON XJONI 
0S°28 00°S2 OS 29 00°0S OS°LE 00°Se OS dl 00° 


9m°O2 CIN) X< (T4N) X) 
01°O-?0HY J TdWuS 

21°O-?OHY JINYL ; 

s2°0:0 : 

Hilbd 3ldWuS ULdAd 


a 


Pee, 





OOT 


0S ‘28 


a secle ei w CU Mestre) 2) eed 


HYIGWNN XSONI 
00°S2 0S°29 00°0S OS “LE 00°S2 oS ‘2] 


8h" O? CCN) X< (1+4+N) X) d 
hl °O-?OHY J dWuS 
hl °O-?O0HY JN 
19°0:0 
HlLod daldWbS blade 


es 





acos 


Ghs par ets. 


69V a fil sauna 


HAYWON X JONI 


O0°SL AIS) Sas. C0. 0s ORS ane UOSSe GS 


6h O? C(N) X< (1+N) XJ) d 
91°O-?O0HY J dws 
hl °O-*OHH JINYL 
EE °O 20 
HLUd daldWbS glad 


6 


00 


224 











ae 


W:0.20d 
TRUE RHO:-O. le 
i. Poe LE eee 0. 10 


-_ SUATIER PLOT 


Z ere eve ee pie | ey 


Ne. *e , : 
i ore e Mg i a 
rare Ie ast pe 
a&y «78 en ° = ry 
a oce ®e ae e 
eS tudv «. 2 oa oe e 
(oO meds. are ‘ 
= Nae yee fe ts Tose ‘ * 
aa pee Cus ge%e 8M : 
Faas tgds fo 28", AiG a 
fe a Re 
He er ice’ sg, ¢ ges A 
d «8 ¢@e 4 Ae é S14 
F ae tee, 
Gt eye ees : ‘te o & 
= phe tier, S "P, .4,32°5% * 
OC Beebe Qh Ge = 
2, ie ey aes iia, ‘ 
-) 


Gu.0 0 sgl gacu u.80 6.40 


BuGOReE, Plt .2.4.10 


220 








|— 

= i 
J 2 
=) 
(Ee a 
a e O_ 

_- 
A ae = 
LE | eo 

(at 
EO) i Ca 
ae 
eo) y= 
Cc) Wt 

= an 
UO rs 
aw ba 
— 
Ci U5 
= Li 
Li | oa 

( 


CO 

= 

= : 

(oO ® 

© ° 

CO Ad 
PE ty \ 





PReUre TLI.E.4.11 


226 


a a 
i a 


BR 4 











"BETA SCATTER PLOT 
a:0.33 
‘TRUE RHO:-0. 14 


SAMPLE RHG:-O.16 
2 EXTREME VALUES NOT PLOTTED 


6.40 





e . cs 
Gin ier ie 
> aa" 7s eee 
aes. 5 
+. ase aay ¢ saa ° 
as " = o ‘ s — 
3 ef t e e e 
a - ~ of eo 
ou? .° OY a, ° 
e $ Pa i ve . 
CO e724 oe Ys, 
Ce weary: ae a 4 
e *e res . eer . 
Ue See “**, ae s e . 
st fais Ne, oe et, : ete a6 
f Ke. sac) we 3 ss A : c e o 
Wide Ee REL al 
+) oe 2 oe cee 
Vibes af Ihe _ ty * 
he rereete ee 
Oo . Pye ee BE ‘ os. 
| et Ae ; ei . * nares : : c 
. = uplie.- ee. Eepel’, JS Sook Say eo as eg te : 
Cc) I 
O00 gaat S20 U.80 6.40 


PLGURE eit E242 


pa ae | 





IV. PRELIMINARY DATA ANALYSIS 


A. INTRODUCTION 

During the period 1955 through 1969 a weather ship sta- 
tioned in the Gulf of Alaska (50°N,145°W) collected, among 
other data, wind speed data at three hour intervals [Ref. 33]. 
The existence of the wind speed data was brought to Professor 
Lewis' attention when a student in Oceanography asked him to 
provide a model suitable for simulating wind velocity data. 
The simulated data was required as input to models of ocean 
temperature mixing. A copy of fifteen years of wind speed 
data was obtained for this thesis. The intent was to doa 
preliminary data analysis and then determine whether any of 
the models discussed in this thesis could provide an adequate 
representation of this data and, hence, a method for gener- 
ating wind velocity sample paths for oceanography simulation. 

Initially, the models discussed here are strong a priori 
candidates for data of this nature. Intuitively, there is a 
strong feeling that an assumption of independence among the 
data is unwarranted. Hence, autoregresSSive and moving aver- 
age models which are designed to reflect the behavior of 
correlated data should be considered likely candidates. The 
non-negative nature of the data mitigates against the use of 
classical time series anlaysis which is based on the assump- 
tion of a normal distribution, and hence negative values, for 


the series. The existence of zeros in the data tends to make 
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the use of transformations, like taking the log, somewhat 
less appealing than otherwise. These considerations indi- 
cate that the models discussed in this thesis should be 
considered as likely candidates for modeling the wind speed 


data. 


B. ANALYSIS OF THE RAW DATA 

There were 43,800 data points available for analysis, 
2920 for each of the fifteen years between 1955 and 1969 
inclusive (the extra data for leap years was discarded). 
Since this size data base made it inconvenient, if not 
impossible, to manipulate by hand, each year's data was 
plotted as a means to promote familiarity with the data. 
The plot of each year's data and the plot of the data averaged 
over all fifteen years (e.g., all data taken at 0300 on l 
January of each year were averaged) are presented in Figures 
IV.B.la through IV.B.lp. Several characteristics can be 
observed from these figures. First, and perhaps most obvious, 
is the expected yearly cycle of the data. Values at the 
beginning and end of the year tend to be higher than those 
in the middle. Second, the data is discretized to a large 
extent. There exist obvious horizontal lines of equal valued 
data. A sort and plot of the entire data set reveals that 
the data consists of values that are integral multiples of 
1.03 with a few values between these multiples. Next, on 


some occasions a series of high values will all be ecual, 
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indicating that some clipping may have occurred (see Figure 
iveremmaevicinity of 2575 and 2625, Figure IV.B.1b vicinity of 
400 and 525, etc.). These last two characteristics indicate 
that statistical properties which are sensitive to the be- 
havior of the "tail" of a distribution may be affected. The 
final observation about the data is that there are apparently 
intervals when the data was not actually collected. These 
instances appear as reasonably long strings of values which 
have a strong linear appearance (as though the values were 
produced by linearly interpolating between two boundary 
values). See Figures IV.B.lh (vicinity 2400), IV.B.173 
Weermnacy 2250), £V.B.1k (vicinity 50 and 1750), and IV.B.1m 
(feinicty 150 and 2725). 

The cyclical nature of the data is somewhat more apparent 
in the plot of the data averaged over the fifteen years (see 
Figure IV.B.lp). Additional evidence of this yearly cycle 
1s presented in Figure IV.B.2. This figure presents twelve 
box plots, one for each month. The data values plotted are 
the monthly average wind speed for each of the fifteen years. 
The interquartile range and extreme values are shown in a 
Standard fashion. As an adjunct to this analysis of the 
year cycle, the coefficient of variation for the monthly 
averages was computed. The coefficient of variation was 
essentially constant. See Table IV.B.1. This will have 
an impact on the choice of the type of model used to model 


this data. 
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Average wind velocity by month for each year. 


Table IV.B.1. 
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This yearly cycle is also shown in the periodogram and the 
log of the periodogram of the data (averaged over 15 years) as 
presented in Figures IV.B.3 and IV.B.4, respectively. The 


periodogram is computed from the data in the following way. 


N 
Let {X,, n= 1,2,...,N} be the raw data and let X = } X. be 
i=l 
the mean of the {X} sequence and as be the variance. Let the 
i, iol 2,...,N; be formed from the {xX} sequence using the 
relation 


where N = 2920 is an even number. The Fourier transform of the 
{Y 3 sequence will have both a real and complex component and 
will have a elements. Let {Z, n = L,2,.0045) be the Fourier 


transform of the rea sequence and let Z., and Z be the real 


jR jt 
and imaginary components of the 5th element of (Zt, respectively. 


h 


Let P. be the 3° element of the periodogram. Then 


ie 2 os N 
: SR -- 257) /2mNo ys J = 1,2,..645 LIV 2S = 2) 
defines the periodogram of the {x} sequence. 

The periodogram dramatically presents the yearly cycle 


(j = 1) as the dominant effect (P, > 150), although there is 


ud 
some indication of a six month cycle (j = 2, P. = 99,0). scome- 
what surprising is the apparent lack of any strong time of day 
effect. The log periodogram reinforces the dominant role of 
the yearly cycle and indicates that six month and six and 


twelve hour cycles (j = 2, 3 = 1460, j = 2920 respectively) 


may be important. 
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The correlation structure of the data is presented in Table 
Iv.B.2. The column indicated as "15 Yr Avg" is the average of 
the values of the fifteen years. The column indicated as "15 
Yr SD" provides the standard deviation of the values about 
their mean. It is not the standard deviation of the average. 
This latter quantity can be obtained by dividing by the square 
root of 15. The last column provides the correlation structure 
of the average data. The estimated correlations remain artifi- 
cially high in this case because averaging reduces the varia- 
bility of the data about the year cycle which intensifies the 
artificial increase in correlation due to the year cycle. The 
correlation structure revealed in Table IV.B.2 for individual 
years closely resembles that of an AR(1l) model, in that the 
k-step correlation is approximately the one-step correlation 
raised to the moe power. The correlations in the table have 
a tendency to be slightly higher than the theoretical, calcu- 
lated value, but the agreement is reasonably good for about 
ten steps. Beyond that point the correlations are kept up by 
the year cycle, which is not as prominent in the yearly data 
as it is in the averaged data. If nothing else the disparity 
between the two correlations is evidence of the existence of 
a trend in the data. 

At €his point sufficient information is available to de- 
termine some characteristics of the general form of the model 
for representing the wind speed data. As noted above, the 


correlation structure is similar to that of an AR(1) process 
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3-HOURLY WIND VELOCITIES; 


DATA: 


The estimated correlations for the averaged data hold up, clearly show 


a long-term trend in the data (in fact the one-year cycle). 


Tableriv.B.2. 
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with positive correlation although nuances may appear as the 
year cycle is removed. Hence, an AR(1) model should be used 
aS a starting point for the construction of the model for 
wind speed. 

In addition, the cyclic nature of the process must be 
modeled. This can be done with either an additive or multi- 
plicative model. An additive model might have a structure as 


follows. 
X = wu. +eEe_, Gly Boo) 


where {Xt is the time series under consideration, LW, is a 
deterministic function of n, and the innovative sequence 

(eee 1s a stationary sequence of random variables. In the 
usual model this stationarity implies that the marginal vari- 
ance aC 1s constant and the correlations only depend on the 


ag (1 2G. , p(X .X o(k)). Using the same definitions 


nee 


the multiplicative model would have the form 
bi Be 4) 


where again the Lene sequence is stationary and independent 
of Wee A characteristic of the additive model is that the 
coefficient of variation is a function of the value of Une 
The multiplicative model produces a constant coefficient of 


Vemslatuen. Since the data, has a coefficient of variation that 
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is essentially constant, in the crude monthly analysis, the 
multiplicative model is preferred. 

We have yet to determine the exact form of the mean, wae 
in equation IV.B.4. However, we do know that this mean will 
have a yearly cyclic nature. We also have yet to determine 
the general structural nature of the innovative process Cie: os 


These subjects are addressed in the following sections. 


C. THE FORM OF THE MEAN; DETRENDING THE DATA 
Two basic models were considered to represent the mean. 


The first was a single harmonic sinusoidal model 


27™n 





= 27mm = 27M 
Pee a ~ by Sin (3959) ~ b, COS (5959) = atk cos (5955 + oe ae 
LEV eG) 
2 272 1, 
where k = (by - bo) and @ = tan “(- 5): The second was 
2 
an exponential sine with one harmonic 
a+b, sin( yh) th,cos (Ut) es kcos (aut+e) 
u = e = ee CE see2) 


The second model has the theoretical advantage that it can 
not be negative and will represent higher harmonics in a com- 
pact form. The sinusoidal model may or may not be negative 
depending on the values for a, Dh) and Do. In spite of the 
theoretical preference for the exponential sine, both models 


were used initially to see if either produced significantly 
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better results. Note that if k is small the models are 
hardly distinguishable. If k is large, the exponential sin 
will clip at low values and is a cycle that would have many 
harmonics in its Fourier transform. 

The values for the constants in equation IV.C.1 were 
determined ina straightforward procedure using the least- 
squares regression procedure of MINITAB and the data averaged 
over 15 years. These estimates could also have been obtained 
from the periodogram at Ww, = ZN; T(w,) a {(b,)* + (b,)*}, 


using the relations 


A N x. _ 
a = ) TT = a iV .C. 3) 
i=l 
s RY Pegs 
oot = 2 ) X. sin(——)/N = imaginary component (IV.C.4) 
. all N 
i=l of periodogram at 
27 Ne: 
a A 27a. 
bo = 2 ) X cos (—~—) /N = real component of 
i=l periodogram at 27/N. 


The varzance of these estimates is 20°/N na iz 9 t= X,'s are 
independent, but since this is clearly not the case here, 
estimates of the variance of the estimates cannot be obtained 
directly. The results of the estimation are contained in 
eouunmn 1 of Table IV.C.1l. 

Similar results were obtained for the constants in IV.C.2 


by a slightly more complicated procedure. In order to use a 
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TABLE . Cd 


Parameter Estimates for Models of the Mean Value 
Function of the Wind Velocity Data 


PARAMETER 


1 hamronic 1 harmonic 

Sine exp sine 
10.236 2.309 
-0.176 -0.011 
22560 0.260 


ESTIMATE 
2 harmonic 
sine 
10.230 
-0.175 
2.566 


=O.o95 


=U Shell 
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2 harmonic 


exp sine 


ENT) 


=O 0 HEIE 


0.260 


=O) 105, 7) 


-0.054 


OU 


- 260 


50S, 7) 


054 


014 


OOH 


AGEL, 





least squares approach, a linear relationship must be estab- 
lished for the mean value of the process. Taking logs is 
the obvious technique to employ, but this introduces a compli- 


cation. Taking logs and expectation of IV.C.2 we have 


E(ln ay) = ike Wea E(ln e) 


j 27n Zin 
at by Sin (4959) + bo cos (5959) +c. 


i! 





For example, if the te) sequence is marginally distributed 
as a unit Gamma variate, G(l,k), then c = w(k), where » (k) 

is the digamma function (derivative of In r(k)). See Cox and 
Lewis [Ref. 29], pages 24-27. The value of the constant c 
will be combined with the constant a in the least squares 
estimation using the ln X'S: giving the constant a' = atc. 
To estimate atc without making Gamma assumptions for the 


innovative process, the X's are divided by 


~~ 


t= : 200 - 27™n 
UW, = By sinls95q) + D2 cos (5959) 





to give Xe The data is then divided by the average of the 


Sam The result of this is a series 


a which estimate e 
with mean value (within statistical fluctuation) of 1 if the 
model for the cycle is correct. The values obtained are listed 


Mee column 2 Of TAble £V.C.]. The results of these estimates 


are in Figure IV.C.1. In this figure the average data 
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is plotted against the value computed for Wn Yommamseth of 
the models under consideration. 

When using a multiplicative model, the residuals are 
formed by dividing the raw data by the mean. The results of 
this procedure are presented in Figures IV.C.2a through IV.C.2p 
using the exponential sine model for the mean. The results 
are not significantly different using the sinusoidal model 
of the means. Hence, only the results for the average data 
are presented for this case in Figure IV.C.3. 

The log periodogram of the average data detrended using 
the sinusoidal model for the mean is shown in Figure IV.C.4. 

A five-step moving average of this log periodogram is pre- 
sented in Figure IV.C.5. The detrending has clearly reduced 
the importance of the yearly cycle, but still shows some evi- 
dence of a six month cycle and six and twelve hour cycles. 
Similar information is provided for the average data detrended 
using the exponential sine model for the mean in Figures IV.C.6 
and IV.C.7. This model does not reduce the effect of the 
yearly cycle as much as the sinusoidal model for the mean, 

» but still shows the six month cycle as being important and 

some evidence of six and twelve hour cycles. 

Since the exponential sine has the theoretical advantage 
of being non-negative and both models of the mean produce 
Similar results when applied to the data, the exponential 
sine is selected as the model of choice and the analysis is 


continued using it exclusively. 


260 








"ME Th. Heb ben? tH wae 





ie 


* 


tHe WI? 


| 
( 


% 


fb HIN? by hd? 





HW A 
ma 


ine tk? 


av ude 








; 


¢ 


_ 





me Ot? 





(> 
2) 


wat Ld 








| /\ 


We Ot 


‘epyep pypustyud GsoG6l 


Uu tel wa ied 








‘ 


UT ITN INOdXx - 


INTM 


ou) youl 


no ust 


Woden xX HOING 


-ez" 


Oh UEt ou Hel 


dA 
IIL 


—S ee 


@ 


“Al 


In} O11 


et 








eanbud 


Oo) Unt 


dheb 00 


ou be 











“ 


/ 


€ 


- 





oo W Qu fy 


nit” 


ho Ut ho oz ow Ul 























| | 


u 


x 


| 





IN TT 1G 
“TPT UU 


uo ue 


3 


42° 


CS 


Gs 


7 
& 


Te 


rd 


9¢ 


pNCiS3s Js3d5 ONIN 


os 2 


x 


in 


te 


“y 


AA yAL 





"pyep popusszod 9667 “GZ° OAL oarnbtd 


stb HAHN & TONG 
wu oud ou nad ol vul wo ued on UHI Hu Ut OU lbat no wet 


‘ 


VU The? no oy? ny ts¢ Ww Uf my Ee mk ee On mld )) het by Ot uth wot it 06 OO ip Oe We ua US OU DS UU Oh UD ae uu Oe ou ot 




















































Ga One in vie oo FS 


3 













vi 
ONIS 


2 
C3305 


we 


Ree 
f ‘ 
New’ 
— 
———— 
wo 
aT 
—— 
ion lo 
a 
—, 
— 
eh fe 
. — 
aa 
es 
—_— 
donee 
‘ 
== 
om 
— 
—_ 


$ 


Nk INTS WEEN ING dX) JAI 
| | 


; | 


~~ 


ot EL tC 


“S 


262 




















































“pyep popuerqzod LG61 "SZ )° Al oinptd 
. 
oe WSHON X ONE 
CULL Sty We tae ime a? on (ae an) OSE Olt Une is Mad a vce uw nle oh One la 08) ou as (ris fhe t ub Oo bu 0 F flu Ob tH HES uo n2t Ou ee ny O04 OD 4h nO Ue DO af G0 ud WU DS DL Dh DO bE OO fi uuu) ou, 





® 
> 
Rs! t 


Gs 


w 


3s 


° 


Si 


322 oc’e 


SIWATISsy G3adS ONIM 


as 


AT 


OS'S 


co 


ages ieee. SINI Sule 


v tes IN J AATIUST Td ELIA 
Cot GILLLOO LIA ON , 


1X | IAW FONTON Wt taf 
7 ee mt \ “UT UCP 


- 


263 








HO vG&c 


i) ud? 


og ure 


ma uae 


mom use 


Du nh2 


Oo Use2 


vies ei eee ene re OEE Oe Cay aa ST a I POE tt MITE eT i TOA fe) fy ren gn WD Bes 1 NE Oe poe. ts at ob ind LA 


"PQeP Ppepucizqod BGbI "PZ-D'AL aanbig 


Ole UAAHIIN KOCINT 
tm O¢e Ou O62 Uf NOE r OT Out ou Oat on ust 00 Oe oo ot oO Oz 00 DN oo Gut i oo 09 oo us i on Q, 


1 6:1 | 8 





a2 G 


ive 


‘2 


- 30 
S’AOIS3e G330dS INI 


$2 


2 


a 


$28 L 


as 


1 SAT eT A) = aNPPRRIIN Balt ait 
| | | ar 


*? 


264 


A we 





"Pyep pepuartzod G6G6T "OZ" O°AL BANbH yA 


(O04 BIOHON MUNG 
LI Dn Get cou 1 oa pet ou O41 on O14 of O41 ou° 














iat din’ att Uae un Gee un old on upd on uRt ou OBl 





Dt Ove Gu O»2 UO Ihe Ou Ned Ou Od 











Ou at Ho ay Ou OS GH OR bi ng ou nz Hout 

































on q, 


~ 
+ 


‘ 


POV TLT HUTA INES (UELNINGdX I JAWLUOE Td EL INN EON TONY ap 
bbtEL SI ttlJollA GNIM J Tk FCT SG TRE e U PUP 


4h 


m 


265 








om Usd ol Gaz 


Ou ee 


uo 092 Ou uu DAL 10 OFZ 


i 


1H} HUTS 
Oo TI 


mw o2e ma Olé 


Hi alecs 


~ 


*Ppjep papus1yd 096[ 


J0o™  WaHKON x10 


Ooo Ou2 on Obl oo vet 00 oit Ou ORI 00 o8f oo Of 








TUT INANOdX I 
JA ONIM 4 


Ni 
co alt 


JATI 
10 





“JC OCA 


00 O2t 


LJ 
| 


OU Olt 


arnbig 


mw Ont on 06 






D1 ld 


-* 


aie 


OU it 


L IW 
| 


a 






ts a oP] 


Oo 


) 





on OS OU O% Ob te 


pa § 


% 








INT CIN § 


C1 | 


lel 


On OL 


nth 4, 


Js°3 


vx 


"2 
S°ONOIS3e J336S ONIPF 


$2 2 39 


Se 


— 5 OG 


os °s 


‘% 


266 








DUK BbRe.* 





uk Owe ud Ute Ltt eh 


( 


\ 


‘ 


‘ 


it le 
Ib tl SAP bLbJorsa 


Ah Ue IMF he Ht Hed 


aie) 


RS Sk 


no teed uu abe 





SIINIES 


“pyep pepterz0d 196T 


“HZ DAL ornbyg 


Oi YIHHEIN WIONI 


UU td ww Oot Ou Opt oo net Oi ugt OO tun! 










on uit 


IU TI NANO dX 4 







eh x ht 


on eel hu o2t On tt mw vot no ar uu Ug 








JATIUOLT Id tt IAW 
HOH & 


TUACETS 






OOS 


Hab Un 












ee 


vu Ot 00 “fhe 





er ul 








6 athe! PRE 1 ate, ty y+ 


uu hy, 


3 


OG * ai: 


aie z a 
Sve001S3Y¥ G330S GNIM 


$2 


38°? 


oss 


ei {2 
i =o 
a, 
een 
a 
— 
 —$e"6 


awa 


267 


+ «Calle 








OO Uac Gud igre’ ut OF2 hi OYE 





OW Hi2 


un O22 wu Hie 


uo hug on ont 


ou o 


‘eqep 


gq! = onan 


PEepus1z2qg ZI6I 


.Obm = =WIBHIN X30NI 
u | 


00 oat 00 i 





Ps | 


"UZ" OAT vanbta 





Ou Nm Gu OCF uw O25 oF ON 08 Ont 


00 us 





uo UF 





ae G 


o 


“ 


i OF wu ul 





G0 h2 















Ou s 


INTUN FHL 1 
“ae 


ts 


8 


38 °S 


oa2 oo as‘t 


= 


2 


00's aL 2 a2 Ss 
3 


2's 


ONOIS2e I530d5 ONIM 


4 


268 








OU Ube (ne ihe 





ut Ore 


i 09¢ WO Uy? 80 OF uo hie 


Ou Ue Ot WIZ 


ANTS 
JU TIA 0 


ta 102 


00 Ont 


oval, al 


OO uel 


N 
N 


-eyep pepua1ied C9O6T 


ott Ord 


al 
| 


On 


uo o4f an of OU Gint 


NO cl 
M 


"TZ°O°ALT oanbhyTy 


YIHHON MINT 


ua “ef oo Of oo U6 












uo art ov vot 


X} FATIUOT Id 1 








oo Un 








oo O04 


TW 


JPL 





it 09 


( 





oa 


\—/ 





ue Os 00 UR 





JIN | 
r 


IN JH AU 
“UOLUE 


cae a’. 
SIMNCIS3e G3s3eS CONIA 


s2°2 


os2 


Qos 


s2°t 


St 


269 





. 4° “| ee 





~ aoe Peter ae 


*pjzep popud1 70d F961 EZ" OD Al ainbtg 


(Ole BIAMIIN XJUNT 
(] ou vAl 


th) O9t 0 oO O84 00 wet oo ot Ou not OU A ou US 


inh UGE im aud mp ate oo uve Ou "4 > One ve DL? oD fed oo ule oo OR mit OBE ou Oi t oo" oe uv §9 





oo ot oo oF i Qt 















oo Opt 














00 ot 
, 












uo 1, 
3 














sé °C 


ea) 


OO Z f° 


> = ; 
ONOIS3Sd Gs30S ONIM 


52 
5 


Py 
> 


oo 6 


$26 


LUOT Td EE AW FON TON THE 1a 
w- & TEDUTS Vee Ulin 


4 


Z7 0 





OW Und ed te bu ore 








| 


| 


a ie 


= "( 


UU thie UU Ort ug es 


higde| | 


~ ~ 


) 


na ace wn Gie 





INIS 
Ji 


"eJep PopuvszVCG GI6T “HZ°D°AT O2anbry 


Ode yas CUNT 
1a Uae au uel OO Eel od ULt oo UT of 056 nu it ed On Uel ou Ob bu dal LO O68 


| 








IST LNINOdXI JATIUST Id I 
JA UNIM 2 INO} & 


AL fie ua a2 OU UY 





1 In 
(1 





W ° 
||, 


Uh us ta Ut La Ut 


INE CIN, 


n 


Ho He ou Ol 








ub Ay 


sf 


ae 


& 


S26 


os't 


Jt I] 
JEU 


OSA 











OO Wier th) Hae 


oo urd 


DO wee 


vo 092 








iw 





Mm 02 


BO O22 









nn ore Ou une Ga 


New 





oO athe OS WAS US eh O84) ee ne .CAINes oF ts te Ee EERO = 


-eavp pepts1z9ed 996T “IZ°O° AL eanbty 


{Ute UIBWAN KONI 
oo DSt OO bet og ngt op n2t 06 OF ao o 








ow Get OW UST 
: t 









Bo Hel 




















a 


INOdX4 SATLUST Id tt 
WH € WATS RE ULE 









_ IAW 


oo 09 











SONTON THT IO 


NO 0s vo q, 


- 
gy C33¢S INIM 


+ 2 = 


38 °O 


ion) 
s7Bndls 


me 4 


os '2 


ass $e°8 ao *s 


oe 


Zz 








OU ted Ou lid im ede 








Lab UNE mh ft 





hu Uke 





Ub ute 





Ou nee OD wie 





JE SIN ES 


Jk 






Ou dud Nt) Wt 





ae 
| /\ 


vu vel 


| 
| 


“pyep popuslqogd LI6T 


{Yl HIOHWHN KIONT 


uo fet on UD) of ol OO weird 








N 
N 


4 


yf) 


RON | oat clr la 
| | 


“uz "Oc AT 


oInbtga 





uFt on et ou Ot 00 vBl 10 ib une fu Of nv oy UD US 


| U 
1H 


JIE LW 


~“ 


4 


2 





W 29 
a Se 


N 


pd yh 10 Oe ou 2 bo Of ita) q, 





GING SENS .g 
al a aR 


: —— 
3345 ONIM 


a8°5 


a 
G 


ao ¢ bi? G&'2 


“52's 


ANTS 





by ved im} Dad ou ii? PUL 





ae 


WS is 


mH) ta im On 


“al BE 


on use 





XX 


09 O22 


Ou We 


ANTS 





-pyep pypuoljyId BI6T 


“UZ TOTAL o4nbtg 


0 HIBGHIN XING 


uu wed Cu Let on o@s on ov no og! 0 Gi 
















On OE 


JU ET NINOIX 4 
TA UNITM 2 


Oo vet on net oO ull Ww (Hp UO fla 00 U8 0 ‘ae ou us 











SI SOE oes 


NOH & TUNA 


o 
2) 


vu os to oh “0 UE Ou tz ata mt ae q, 





ON TON REL HO 
UT UE 


sa 


sd 03305 


3 


O56 


i 


Oni 


OG 2 


322 
Siehais 


“2 


ws 


e- 


274 





ih Ube 





mb tied 


m O12 na OY? 











be 


mo 


| 


00 Uce On Ul? 





OU Ge 





TA INIS 
JU 





uo un? 











"PQep pepttv1790qd 696T 


Ob ORHAN K GONE 
oO Vat ow O71 Gu B91 ob ust 0 0 it ee | 








on ust 








TULININOdX I JAI 
UINTM A hit 


ou Udl 





| 
| 





| 
| 


oO tui 


‘OZ° OAL Banbia 


On Ubi 0 06 On OW OO ne un o9 





















JO) We TERA 


nes POT La & 


ou Os on & 





Ho he 





Ou n¢ 























os: 


bas 


‘2 


20 


“2 


ONDiSSe GiadS ONIM 


fad “6 


27> 








And dtodd 





15 





tind om hid eo iw ULé 








Pa 
Ue TAU 


uu Une 





alla. 
SLL ape 


AN 


s 


0 G2 OO Gee wo Ute 


fe 
| 


(ou US 


OO nde ca 


») 


| 
om 


U 


| 


Le EyUBUOUXD) 


| 
I 


Oa en 


popu yop YPea vbULaIAY 


st YBMEIN ¥40NI 
vd O8! 80 O48 8 86—DGl te 








ia) Gat 


a2 Sl 


oa yl oo a6 


aaubt 














(INTIM J TUNOH-& 


oo at vu°os oy As on On OO bE 





ININOdX4 SALLVBOL Td EL INW FONTON I 
| 


TC Sel et 


ou a2 





tur od 





Gina 


"2 


a2 


39 


ue 


WWNCL3s3e J35dS UNIM 


- 
> 


as 


276 








fat Ue? uu owe Mm He om used 





be 
MU 


| 


i Uae? 


| 
/\ 


| 


| 
J 


ou uv 


ie 


on Le HO U¢2 mb Ole 


aoa ee 


het hI 


to Oe 


As 


a ee a _ 


*(auTS) pepuarqep qed aHeABAY 


I] 
orn | 


tht OBE 


l) 


no Owl 


ow o1t 









SUNTISS 
(1IO VA U 


uO opt 


sot  wIHWAN x IUNT 
bu bt 


00 OG t 





dA 


L | 
N | 


00 O&1 


US 
M 


"E*O'AL oanbtg 


o0 Het 





} | 
I 


| 
| 


00 olf = ow vO 00°06 on 08 


dil IMW 
5 Diets 


Ho Ot on ny 00 as no Of Ou OF ae uo ul 


SIN Nels iG 


WAI SAY -ULUC 


ou 1, 


is 


$f 


G336S ON 


ua 


me 


“f 


ams 


“’ 


ii 








| J 
ere tele 


*Ssodoad VATSSethasGyNe AVps1o-ysary SHLTQuasoT wNAods punosr by AeA ‘apodAn Aparvod jo FRAoueys1 
49})}e SBTOAD ANOY BATEM] pueR 


el A We 4 
cll}. | o 


XIS JO BoOUSSOId puR apPOAD YYUOW-XTS JO SouvuTMOp SMOYS ueTHOpo! 130d boy oO OA 


se a 16 Hud é : Ih a ut 


-e JUNWTS ye 


# 


INES JATLUOT Td TL TAW -ONEON Felt ot 
| 


with 


og Q 


ee) et 


oc *a- 


30° 


Sale! 


a5 


pee 
SVIBA ABHIDOSIs 


+26 


oo 


00'S 32 


30°% 


~ 


ac 


23 








: 


| 


| 


WI NU 


HOURLY 


w2 


Ja°ie 


80825 30 
Sn SA wm 


2 


or. 


su900 


79 


o3°5° 


00° 8- 


oo°e- 


GO°e- 


a2 2 he 


i] 


a. 





# 99 


4 0 Ay 95 06 16 1 2? 308 1.49 
PRE QUENT) 


64 


n 52 


aa 


QO ua 


. £ 
ace ve 


Figure ivV.¢ 


a er 








IIs 
Mel JAY 


sazpodAo Apaevd jo TRraowosas 
reqzge sapoAyd anoy sapomy puvw xXts JO DOUaBOId pue OTOAD WYUON-XPsS JO GoUeUTMOP SAMOS weabopotszod boy 


EN Wee Teh 
6s i 6h i ae i i2i vi wo 4 ob 0 he Oo té &3 0 


a me re ee 


| 
‘| 


' 
| 


’ 


I 
i 


+ 
{ 


i ’ 
! 
; q 


OTS Saal NS ey 
ae eal | 


———= — 
/ X 
ee 
N ‘ 
\ ome — 


JA GNIM JL UdNOH-& TUR 


N 
WU bel tel Ob 


“Sa ont 


DdIXd AATLUDT Id PL INW *ONTON Tt 


ale 


oanht4 


| 
=m) 


——_e 0 SS 


(| 
Lf 


30 


33° &- 


33°S- 


be T+ e- 
aN oA wmowIO0Slyje 


a7 


ac 


“Je 


5 
d 


vaz"i- 


oe 


4 


a 


280 











J 


| 





} 


| 


XPONENT TF 


J 


HOURLY WIND 


Bi 


. 
V+ 


» LUA 


PRE NU EN 


ine 


bee 


Moe 


yt NI 


( 


a 
o 


) 


SSS ee 


90° 2- 
SNA 


CC °E- 
woe5O00 


261 


che 


ledge 


__—__————— 


20°3- 


30°S3- 


- 
= 
- 
fal 


13 


? 


2 






sa 
PHEQUENEY 


IV 47. 


ie 


o 20 


0 04 


>. 


* 


Figure 


D. RESIDUAL PROCESS PROBABILITY STRUCTURE 
Having removed the dominant seasonal effect from the 
data, it 1s possible to investigate the structure of the 


residual process =a in the model 


X = UE€_. CiVveD . } 


The two facets of the probability structure of the stationary 
process te, which were addressed in Chapter II were the 
marginal distribution and the correlation structure. The 
residuals produced by dividing the raw data by the appropriate 
value of the mean were supplied to HISTF, a histogram and box 
plot routine developed at the Naval Postgraduate School. 
Histograms for each year and the entire data set were pro- 
duced. These histograms are presented in Figures IV.D.la 
through IV.D.lp. The shape of the histograms is consistent 
over the years and indicates that a Gamma distribution is 
appropriate for modeling the innovative factors. The param- 
eter k can be estimated as the reciprocal of the coefficient 
> of variation squared (see equation II.B.4.16). The estimated 
value of k for each year is given in Table IV.D.1. 

A careful examination of the statistics associated with 
the histogram will reveal that the values for the skewness 
and kurtosis are low compared to the theoretical values for 
the Gamma distribution, namely 2//’k and 6/k, respectively. 


However, this is not unexpected when one recalls the 
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Table IV.D.1. 


Moment Estimate of the Gamma Shape Parameter by Year 


neaw Estimate of k 
LAS Bess) 3.96 
ILS Se 4.16 
1957 4.38 
1958 3.68 
Les ey, Sale. 
1960 4.45 
1961 cleo, 
1962 Seo 
sie 3 32,86 
1964 4.49 
1965 3 Z 
1966 4.04 
ood Ss) glee. 
1968 4.24 
Ison, Sra: 
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discretization and clipping that has apparently occurred 

while the data was collected and processed. As a rough 

check on the extent of the clipping required to produce 

values for the kurtosis and skewness similar to those for the 
data the following procedure was examined. A sample of 2920 
values for a Gamma distribution were produced with mean 1 and 
shape parameter k = 4.0 by a call to the NPS random number 
generator LLRANDOMII(SGAMA). A histogram of these values was 
produced sequentially for the following cases. The data was 
clipped so that all values over four were set equal to four, 
all values over three were set equal to three, and finally the 
highest ten percent of the data was set equal to the value of 
the 2890 sample order statistic. The first four central moments 
were estimated under each of the conditions. The results are 
presented in Table IV.D.2. The results indicate that a 
clipping of the top ten percent of the data will vield results 


for skewness and kurtosis comparable with those observed in 


the data. 
ABIDE LV 6D 
Mean SD CV Skewness KUreEOSsis 
Gamma 0.986 0.504 aim Sale eS 2 lea 
Gute at 4.0 0.986 0.504 Ore Sekek eee2 3 Zens 
Gime at 3.0 0.984 0.498 02506 0.994 IesAc6.5 
bo. Cut 0 Sic 0.489 0.498 0.861 O2516 
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The conclusion that the innovative factors can be modeled 
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Since the estimated correlations o0(k) are affected by 
remaining trend (as seen in Table IV.D.3), it 1s best to 
examine the structure of the dependency process via the 
periodogram. 

Figures IV.D.2 through IV.D.5 show the periodogram and 
log periodogram for the 1955 and 1969 data detrended by the 
single, yearly harmonic exponential sine (see equation IV.C.2). 
Superimposed over these plots is the spectral density and log 
spectral density of a theoretical AR(1) process with correla- 
tion equal to 0.849 (see Table IV.D.3). This spectral density 
is also the spectral density of the GLAR(1) process. We have 


2 


eae ae cos (w) 


with »9 = 0.849. 

All of these plots show that the detrending has reduced 
the importance of the yearly cycle and that a six month cycle 
has now become the dominant factor. The theoretical GLAR(1) 
spectral density fits well for the periodogram after the 


point representing the six month cycle. The six month cycle 
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High values for average data shows that a remaining trend is inflating these correlations (see Figure 


TV2D27 2) 


Table IV.D.3. 


1 HARMONIC EXPONENTIAL SINE 


DETRENDING: 


3-HOURLY WIND VELOCITIES; 


DATA: 


a er es ee 


15 YR SO AVG OATA 


15 YR AVG 


LAG 1955 1956 1957 1958 1959 1960 19€1 1962 1963 1964 15965 1966 1967 1968 1969 


ee Se =< > —_— See a — ee —_ = <——- eS | — — —_— = 


——— 


PANG PFPPOT SUVA QOD PHP ONOM ANWNAHOWM SAITO 
FPTDODON HOF IHD DWI HN OO OIA AAI SF SW Pe 
DP UV SP PPO OIC THAIN OY QUIEN QUAERN QUIEN ON AION CQUIGQUINGQUOQIEINY 
**eoeeeeteeeee © 6 & @o © 8 oF 8 8 8& @ 8 6%  G 


YVOOODO9D999D9 909900000 CO009D0900000 


O POD DQ O44 Od dh QP mt Bet HOA OE SOUMN HOO CICA S 
NAAN AI SF PPP PP TOPOL PPP PEM QM MN 
ODO9900 090090 00000070 090900900900 


eoeere ff fe Fe Ft @ © eeeeset © 8&8 © © © 8 & 8 6 8 8 


TFQ9OT9D099000009090090000900000030000 


QUAD PFO DD MAE FD ANIDM MGM OW 00 TWAS BHAA Pe 
NEN & DAM TAMA SCN NNO ODOR DDMOM Qonunmwnd 
D OWS MOON Nat att at St et IS HOO OOD OOOOCOO00O0O 
*e © 88 tt 8 © 8 8 tH HGe ef Fee eH Feeee sees 8 @ 


ODAQOO OSS QW OO OOS OSS OVO CSC COOCOS000 


WW HDD MD aO MO OD DAD pre 9 Ft OD UN AICI RCI CSP 
MW CAM AICIN Nae 4 tS Ht Hat HSI HOO ODOC OQOOOOO 


*e* ¢©¢ ee? © @¢¢ ee © Oe © HF 82 Ff 8 Ff 8 6 8 6 8 GS 


DOO0C0O900 90 000000089000 CO00090000000 


DOLPOAPOOM LFHAO HOM DOD DHMH PHO FONORDDNOA YW 
BOM QUE MAI OLIN ad tet et hed et ed at Ht et Ht Ht HOO ded ed 
®*esv8te8 © © fF 8 6&9 8&8 8&8 8 8 oe &€ Fo & Fe BO FT Oe Fe 8 


elolelalelelelelalelalealelalelelalelalalelalol*lalelelelelalel=) 


NIDDM ALFRDO DMO QM Wee PO PE Om Or FOI RY 
DOWS MMNUAIN HS HOOOCODOOODOCOCOOCOOCOOO 


®eteertkeesetee8e8eee 8@ 888? &©e © © © eo 8&8 86 @ & 8 8 


OOQQMOOC OO CQSCO9OOOOS SG CCO C9 0C000O90000 


PODS ADD DOD PF MODD OM MANAINN NCA MH 
DMOMP OM NA aa tt OOT ODODOOCOSGOCO09D00 


ees oesetete © @ @ @@ &©&@ees#t6t & © F&F ee ec OhUhTlC HhUchHC HhUhOlhUhHmhUh}Hmh(CU }HlhUF 


OO90000 009000009000 000 O900090000000O 


DD ARROAIPS CAD Pi FNAL SHO HIN COCO COP LP OCIAQI RHO OND 
BRINE PONV Ate tte et tet ed aed et tt mt ent ad ad OC 
*eewteteettses2seses#ee?t#eeete#e@8es9 @@ © @ @ © @ &@ @ @ 8 


STOO OOTOO OOS SOOT COO OCO S90 0COOO9D0009 


fom ors dvatatoolonl@si atoalacl?ola tol Sin ol ahs ae maakt oom dalla fevletr rola! 
DOW-F OANA AHO ODO ODO OOOCOCOOCOOCO009300 
eoeteoeeeseecewvrs#s##eeset eet eFeertkeeetee* @ es © e © ® © & 6 @ @ 


DODO TCOQVOQV OO OSDO COO OOO COO90COS0000 


PHAM POM TOD DMNA LN Halt HHO MON =-1 ODM OO 
DPM WAI CIR let at ett et tet St tad Heed HI SH HOMO O 
oot e*ee ee eee fe Fe %Fee ee © Fe FH © F&F © 


DOOD TO OOO OO CTOQ 0 OT COOH SSOCCSLe0ees oY 


ANID DHA SAS Poh SO Pe Pe Pope DDI PPM Rt Bh eu or PCO 
DMAP OMAN NA Ae HH HHH I HOODOO OTOCOOOO 


eeeeseeeeest?e 8 @ © ee © 8 @ & & & F&F Bef F&F 68 & @ 


O OVO OOS OOOO OS OOO OO OSS TO OSSOCco9d0000 


NE DOUALA 0 WD DOIN OF FOO SN N et Nat Se NI Sat 
DOWIE MAN at OO STOO OOOCOC COSCO OCOCoUoo 
*e*8ete @ ee @ ee ®t 8 fe Ff GH eee fF 8 6 Fe Fe oe HB & 8 


OWDSOCOSOSCOSC COPMOO TOC OTC VTOCeCcooo 


OM QO MOM AMM DOS MAIN AN COCOA AIS PLAIN 
DOWPFPMOMNAIAINISHHaOOCCOOOOCOooOoCCcCoscceceo 


ee2oeoeestse8s @ ge &@ ee @ &et¢ Gee ec emhUcCOhUCUC MCU OFhUcMhUhO FH FH HH SF HS 


DODO TOVOOCSO STO STOO CSO COC eocoeooceoo 


DHROOFRHROMODOM OOD OF FLOINN He CHROCORW 
THM ALAPMOAICIN NA tet at Het et Het Het Ht tt HD et tt ot ed 
*es28k8eegees#s8ef tees fF 8 c8OmhUMHMhUhO TF HH FH FH HeHee FE Fee 8 


oo tmooccoccoocece oc ooeoceococroooeo 


NADONMS HDF NAHOOODM™ OU LROAMMNOOOCmaM OM 
DM OWNS COC MAN AIQUAICIQIN Qld at Ht et tet at eg Ht HI COCO 
® @©@eeteeefiet te? &@h—Uh ohUC<C DCM WC }OhUch OHOChUC HmhUCUC }H HUM HOU HMmhCUCMPMMUCUCUMhOChUCUCMOmhUCUCOhC<CHMMhC<Cwh hCUhThUCU MHhUCUhTChClUh HhUCUcFhlhlUhH 


DVOOOCOOCSCCODWOOOCOTOS COOSSSCECWEOOOO 


=P PFW OCD FOP HO DAMS CUICS SIMO ANI O ODD OWING 
OD OUD SP CAIN OA St et et et ted et HH HR IMHO OQ OC OO 
eet egveeeeetrteee tt ©? eteteteeeee eee et eee ® 


OTOOCOOCOOOCOOSCTOSCCOOCSC OOOCOCOCOOoO 


OOTP PFOdDWN NY DIN OIA OTA MOANMAIMN AICO @O@Oruw 
DULY SFP VT NS Nt tt th tt et tet et et ot tt St et St It CO OOO 
eeeee8eet® eee? @e 8 Fe © 8 F&F © Ce FF © Se B 


oooocccocoovoecEeceoccooqoecoce ccoceo 


CIP OWNO DAI DLP AIKQaMOnANOA OM OO OWOOINAS — Fn 
DOP MAUNA Hat Hatta ODSOCcecooeececcoo 


oe OP Peesrsetseee#rnfseeerke?eseeseesete © @*# @ 8 8 © @ & 6 @ 


qoo0o0ccocococcooooooec ecocoocooecccoo 


HNIC FINO DRHOANOLAOE DRO AINOATHAOP DTO1 
tet tet at tt eH HA CUE AION OI AIO IONE C9 6 


602 


eee a ee a ch A CL ES A cee tl A ea a IE LS I, AS 


— — a i- — a = 


— a a eo SoS ee ot — — ome 


= — a = aes ee 


SS OF EE ER a 8 em eee ~ 


Sed 





eu ¢ 


alt 


ouowIey T Jo wed 


Jemiciiunete! 


mor 


pee 


SelGeie a) 


1) 


2 


ceé 


“6 I 


AH OTNCE DY 
es 6h 


boportad 1aA0 pesodwitedns st 678°O JO WOLTeyT 


4 Lis 


t 


Bi 


ont 66 0 


+4 eR 


rir were Ty) 


9ii02 YAM ssooo Id 


aida fr 


bal 2d yep popuol VYVp 
(TL) uv 


he og 


yo wnayoads “7° 


Ly 


ave Lise Sadi Ua dU WM i Whi i i} fi 
i 


ao ty Oo 


ee 0 th O 


° At 


th 0 


oanbtd 


bn 0 


<n C 


SINT ONS =t0) 
( 





€a 


8 


oc'2 


0d *& 


AS 


ie) 


“a 5G 
9003tus¢d 


a9 


*§ 
SNA wee 


GGkL aoe: ao*ts 3o°ol oc 


s 


ao°h 


si 


. 


oa 


oo°sL 


303 





—— 
rd 5 
pect’ 
ion j 
Jen zoey 
a 7 
+, { —_7 
ce 
~— 2 
=e 
i= | 
Si oN 
poe ee 
i 
=e 
(at 
F teal 
SS 
—— 
f= 
aan oes 
oo 
(Va 
. 0 Meee 
= 
——s 
——* 
———————— 
oo} a's woz 


—' 
i 
fom 
: a 
— 
_~ 
New , 
— 
— = 
fy ee 
—a 
1 ——! 
tT —e 
oa oe 
=—— 
alee 
7 
a 
— a 


7 x 
——/- atmo 
K fas 
\aee/ 4 
re 


ee ee a ae 


304 


o0°8- 20° n= 30'S- co 
SNIwbA wedDOOO dade 3907 






A 


A 


Go° B= 


; 7 


OG Qi- 


= = = ® i 7 “ “ 
oa Tee: - cs re Se 
ES. Eee ale faerie is wath 4 


id Magri pb eG TB 
MO nek 4 std gO Aa 
- =a » a = * : 





ay 


2 We 


f Oo? 


54 


FRE QUENET 


Figure IV.pD.3. 
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has now become the dominant factor. Lacking in these plots is 
an indication of a time of day effect. The appearance of a 
time of day effect is limited to those plots which use the 
average data (see IV.C.4 through IV.C.7). Figures IV.Cc.4 
Mmomeugm LV 7C./ and I¥V.D.2 thorugh IV.D.7 indicate that a 
further refinement in the model of the mean to include a six 
month cycle may be helpful. This topic is covered in the 

next section. 

The correlation structure of the detrended data is depicted 
Pm er lveprs and Fagures I[V.D.6a through IV.D.6p and IV.D.7. 
Since the one-harmonic year cycle in the data has been reduced, 
the correlation of the average data in Table IV.D.3 more closely 
reflects that of the average of the fifteen yearly correlograms. 
The higher values for the correlations in the average data and 
its failure to fall below 0.20 may be indications that a trend 
still exists in the data (the six month trend) which is arti- 
ficially inflating these values. This may be a further indi- 
Gar ioOumou cme Gesirability of including further cycles in 
the model of the mean. The slight increases in the correlations 
for lags of 8, 16, and 24 in the average correlogram, Figure 
Iv.D.6p, and the correlogram for the average data, Figure 


IV.D.7, may indicate a small time of day effect. 


E. REFINING THE FORM OF THE MEAN; A FURTHER DETRENDING 
Since several plots in the previous section indicated that a 
six month cycle had become the dominant factor since the removal 


of the cne-harmonic year cycle, a further refinement for the model 
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of the mean was considered. In this refinement, the sinu- 
soidal model for the mean was reintroduced to the analysis 
to see if the addition of the second cycle would allow this 
model to outperform the exponential sine. With the two 


cycles included the sinusoidal model becomes 


Z : 27™n 27n 27m 
~ a 
A a + b) sin (3959) + by lsg55) + 63 sin (Tyg) 
+ b cos (452) CTV |. 
4 1460) ° -E.1) 


The exponential sine becomes 


27n 27n 
atb, sin (33%) +b, (sayy) +b, sin (F7gh) +b, cos (T7¢q) ° 





(EV... 2) 


Parameters for these models were determined following the 
procedures in IV.C and the estimated values are listed in 
Table IV.C.1. The plots of the two resulting values for the 
mean are presented in Figure IV.E.1. Since the two curves 
_are nearly identical and the exponential sine is preferred 
on a theoretical basis, the sinusoidal model was not further 
considered. 

Figures IV.E.2 through IV.E.5 show the periodogram and 
log periodogram for 1955 and 1969 after detrending with the 
two harmonic exponential sine. The value for an AR(1l) proc- 


ess is superimposed as before with a correlation of 0.794 
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(see equation IV.D.2). All of these plots show the yearly 
and six month cycles much reduced in importance. They also 
show some weak time-of-day effects, but these are more 
noticeable in the 1955 data. 

The correlation structure of the data is shown in Table 
IV.E.l1. The average data correlations are now lower than 
those of the average of the fifteen yearly correlations. 

They also drop more quickly than that of the average of the 
fifteen yearly correlations and eventually go below zero. 
This is another indication that the trend in the average data 
has been largely removed. Figures IV.E.6 and IV.E.7 show the 
periodogram and log periodogram for the averaged data. As 
has been noted previously, the time of day effects are more 
noticeable in the averaged data than they are in the data 

for a single year. However, the effects are prominent enough 
to warrent further consideration. This subject will be ad- 


dressed in Section IV.G. 


F. RESIDUAL ANALYSIS 
Since a first-order autoregressive model appears to be a 
good fit to the innovation process {fe}, we need to examine 


this hypothesis more clearly. If we were dealing with a 


linear AR(l1) model for the residual process 


= = p¢e ate. Clan Jy 


where Yn is a sequence of iid random variables, then computing 
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as an estimate of pe. would be of interest. The estimated 
et Should have a flat spectrum. Using the Gamma generation 
scheme of equation II.A.5 (i.e. = 

q cnc ce, fn = Baten-1 + C,S,) reduces 
the value of differencing since the coefficients in the genera- 
tion scheme, Ba and Ch are continuous random variables and 


not constants. However, this differencing procedure may pro- 


duce some insight to the data. Hence, the differences 


x one o(L)e, 4 (Tie es) 
were produced, where 0 is a one-step (lag one) correlation 
for the two harmonic detrended data and En and En-] 2re two 
harmonic detrended data values. 

Since the data has been detrended and the differencing 
serves to remove the dependence from the data, one would expect 
the periodogram of the detrended, difference data to be flat. 
The periodogram and log periodogram for the detrended, dif- 
- ferenced data are presented in Figures IV.F.1 and IV.F.2. 
With the exception of a relatively strong indication of a 
Six and twelve hour cycle, the periodogram is, in fact, 
reasonably flat. The log periodogram indicates the same 
characteristics. 

The correlogram for the detrended, differenced data is Figure 
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are relatively low, indicating that the dependence structure 
has been largely removed. Second, the alternation of the 
slemmeot the correlations is an indication that there st#1l 
eXists an important cyclic component in the data, in particu- 
lar an alternation of twelve, or six hours. Differencing 

two sine functions (i.e., sin (422) = sin (=2ntt)) will pro- 
duce a cycle with period of two if they have non-zero ampli- 
tude. Therefore, the alternation of the correlations is 


evidence that an important cycle still remains in the data. 


G. A FURTHER REFINEMENT OF THE MEAN; THE LAST DETRENDING 
Since the evidence of the preceding two sections suggests 
that there is still one important cycle in the data, a further 
refinement of the model for the mean was undertaken. The 
evidence suggests that there may be six and twelve hour cycles 
in the data. These cycles may be the result of the passage 
of pressure fronts over the data collection location. 
Only the exponential sine model for the mean was con- 


sidered in the final detrending. The final model for the mean 














was 
27Nn ? 27mm 
i. EXP [a+b 18in (sag) + byc0s (Fggq) + by sin( 4) +b,cos (S775) 
a “ts — 
eV ae eal 
4 6 } Gin Ged) 


One should note that the sine function with a period of two 


is omitted from the model. This is because the sin (nT) is 
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identically zero if n is an integer. The implication of this 
is that we essentially lose the ability to determine the 
phase shift for the cycle with period two. This may mean 
that our attempt to remove the six hour cycle will not be 
completely successful. Parameters for the model in IV.G.1 
were produced by the same techniques used previously (see 
Seeeron £V.C and Table IV.C.1). 

Figures IV.G.1 and IV.G.2 show the periodogram and log 
periodogram for 1955 data after detrending using the model 
of the mean in IV.G.1 (see also Table IV.C.1). With the 
exception of the six hour cycle, the periodogram compares 
favorably with the theoretical AR(1l) periodogram superimposed 
Over it. The log periodogram shows the same characteristics. 
Similar information is presented for 1969 in Figures IV.G.3 
and IV.G.4. The strength of the six hour cycle is reduced 
for this year. Finally, the periodogram and log periodogram 
for the averaged data are presented in Figures IV.G.5 and 
IV.G.6. The comparison of the averaged data with a theoreti- 
cal AR(1l) process is considered acceptable. 

Note too that in Table IV.G.1 the 15 year average 
correlation is commensurate with the correlation computed 
from the averaged data. Thus the discrepancy between these 
quantities noted in Table IV.B.2 has been explained. 

temas, Se worth noting in passing that a SUrprisingmsresult 
of this analysis is the failure to detect any multiple-day 


cycles. Apparently some meteorologists believe that there 
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tHe QUE NTE F 


IV.G.6. 


Figure 








1s a six-day weather-cycle driven by the passage of storms. 
This analysis has failed to detect any such cycle. It may 
be that the high correlation among the data and the expecta- 
tion that the actual storm cycle will be reflected in the 
data has created an impression that these cycles exist in 
Piemcaata when, in fact, they do not. This contusion of 
quasi-cycles produced by high positive correlation and com- 
pletely deterministic cycles is common in applied science. 
Figure IV.G.7 shows a sample path for a GLAR(1) process with 
high correlation, p(1) = 0.85. Although it may be tempting 
to conclude that this process is showing evidence of a 
cyclic nature, there is no deterministic cycle in the data 
shown. The behavior displayed in this figure is typical of 
an autoregressive process with high correlation, and no cycle. 
Apaplemseeacorrelations for the 4 harmonic detrended 
paeoeceesevlded in Table IV.G.1. Its characteristics are 


much the same as those of the two harmonic detrended data. 


H. SUMMARY 

The model suggested for the representation of the wind 
speed data now has the following form. The basic structure 
is that of a multiplicative model, that is it has the 


form 


X a Wet ere oe ie rr ore c (IV oH) 


347 


—_— 


- 
7 





UOTIETSAIIOD YHtyH ULTIMA ssac0ig 


SATSSOIbsrOINYW ue jo SANZeEN ,OTTOAD tTsenG, SMOoYS YAed oTdwes L°D °ALI FYNOIA 


: YIGWNN XJONI 
oo°00e2 oo°s¢ét o0°0ST oo'’Set o0'DD! 00°SL 00°0S 00S 00° 


0. 
CS 
an) 


pst 
iene Lica ay 


Oe 
IN 1d 


S@°O #OHY JNYL ds 
9°0=0 'O0°h=)y 
HluUd 4aldwWbS (1) db 19 


O39 


7A 


OS 


348 





Table IV.G.L 
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The DS sequence represents the raw wind speed data. The 


u. is a deterministic function of n. The innovative terms 


n 
ad are modeled by a GLAR(1) process. 
The GLAR(1) process was discussed and analyzed in Section 


Ii.B. The generation scheme presented in equation II.B.1.1 


1s repeated here (with oa replacing X.). 


= = A € se BAe : CLV) 


The innovative sequence {te} is itself correlated. The 
parameters of the GLAR(1) process control the correlation 
structure of the model. (See Section II.B.2, in particular 
equation II.Bl2.1.) 

The mean es has been modeled as a four harmonic exponen- 


tial sine function, 














a .; 27n Papel ; 27n 27™n 
Se EXP [atb, sin (5955) +b,cos (5959) th3sin (7GEq) tb, cos (F7Ep) 
27™n 27m 27 
+b, sin ( 7 ) +b cos (—j7—) tb.cos (=5-) | 


(Cli soy) 


The four harmonics represented are a yearly cycle (coeffi- 


cients b. and bo), a six month cycle (coefficients Db. and 


1 
Dads a twelve hour cycle (coefficients De and be) and a six 
hour cycle (coefficient bo). The values for these parameters 
and the "a" parameter are found in Table IV.c.l. 


The innovative terms are modeled by the GLAR(1) process. 


The parameter values for k and q were determined to be 2.843 


S)3)10 





& 


and 0.727, respectively, by using the numerical approximation 
to the maximum likelihood method described in Section II.B.4. 
The data used for this evaluation were the residuals produced 
by the single harmonic exponential sine model of the mean. 
The parameters were not recomputed for the two or four har- 
monic exponential sine model because of time limitations. 
These parameter values give a correlation of 0.744. This 
1s somewhat less than the average correlation of 0.826 EOE 
the single harmonic residuals (see Table IV.E.1). However, 
this deviation is not considered serious. This is because 
the estimates produced by the four harmonic detrended data 
may differ from those produced by the two harmonic data and 
the correlations for the one harmonic data are modified by 
the presence of the six month, twelve hour and six hour cycles. 
The simulation study of Section II.E indicated that for 
large k and high correlation the standard deviation of the 
maximum likelihood estimates was about half that for the 
moment estimates (see Figure II.E.1 and Table II.E.2). In 
addition, neither estimation procedure had any apparent bias. 
For these reasons the maximum likelihood estimates are pre- 
ferred over the moment estimates in this case unless com- 


puter time is limited. 
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